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to  40  km  range.  Events  listed  with  hypocenters  at 
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depth  Interval,  with  the  stress  levels  being,  again, 
near  the  1 kllobar  level. 
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at  the  shallower  depths.  It  is  likely  to  be  the  depth 
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Eplcentral  stress  map  from  events  in  the  Japan-Kur iles- 
Kamchatka  region  for  the  depth  range  from  20  to  33  km. 

High  stress  zones  occur  along  the  Kamchatka  and 
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of  zones  of  recent  large  earthquakes  such  as  event  no.  61. 

Eplcentral  stress  map  from  events  in  the  Japan-Kuriles- 
Kamchatka  region  for  the  depth  range  form  33  to  40  km. 

The  high  stress  zone  off  the  Honshu  coast  appears  more 
regular  and  extensive  in  this  depth  range  and  constitutes 
a high  risk  zone — especially  in  view  of  its  proximity 
to  Tokyo. 

Eplcentral  stress  map  for  the  Japan-Kuriles  seismic  region 
from  events  having  focal  depths  in  the  depth  range  40  to 
50  km,  with  events  occurring  before  the  large  event.  No. 
193,  represented.  This  event  and  its  aftershocks  occurred 
within  the  moderately  high  stress  zone  off  the  northeast 
point  of  Hokkaido  and  results  in  a considerable  change 
In  the  stress  level  and  spatial  pattern. 

Eplcentral  stress  map  for  the  Japan-Kuriles-Kamchatka 
seismic  region  from  events  having  focal  depths  given 
as  33  km,  with  only  the  events  before  the  large  event 
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Figure  35.  Eplcentral  stress  map  for  the  Japan-Kurlles-Kamchatka 
from  events  with  focal  depths  at  33  km,  with  event  No. 
144  and  Its  aftershocks  Included.  The  stress  field  in 
the  northern  high  stress  zone  Is  now  much  reduced,  with 
some  high  stress  Indicating  events  occurring  near  what 
appears  to  be  the  edges  of  the  new  failure  zone.  The 
stress  level,  as  evidenced  by  the  event  sampling,  there- 
fore appears  to  be  much  lower,  as  would  be  expected,  but 
with  higher  stress  near  the  edges  of  the  recent  failure 
"plane"  formed  by  the  large  event,  as  also  would  be 
required . 


I,  Introduction. 
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The  principal  objectives  of  this  study  were  to: 

(1)  Refine  and  extend  both  the  theory  of  earthquake  radiation  field 
generation  and  the  approximate  models  derived  from  this  general  theory. 

(2)  Apply  the  theory  in  a manner  designed  to  provide  a iheoretiial 
basis  for  the  currently  used  discrimination  methods,  in  particular  the 
“b  Z®.  method. 

(3)  To  design  and  verify  new  discrimination  methods  which  arc  simpller, 
more  reliable  and  that  apply  over  a larger  magnitude  range  than  the  ^ 
method . 


(4)  To  reliably  estimate  earthquake  source  parameters,  in  particular 
tectonic  stress  drops,  using  the  theory  as  a framework  for  the  interpretation 
and  as  a basis  for  defining  appropriate  observational  measurements. 

(5)  To  study  earthquakes  found  to  be  high  stress  drop  events  with 
emphasis  on  testing  the  discrimination  methods  with  these  events. 

The  results  of  the  work  in  these  areas  is  discussed  in  the  following 
sections.  The  emphasis  has  been  nearly  equally  divided  in  the  first  four 
areas,  while  the  final  objective  has  been  much  less  developed  in  terms  of 
final  results.  The  theoretical  work,  discussed  in  the  next  section  11  has 
resulted  in  the  development  of  a very  general  theoretical  framework  for  the 
description  of  the  radiation  fields  from  earthquakes.  We  are  still  developing 
specific  models  involving  considerations  of  failure  zone  geometry,  boundary 
condition  approximations  and  prestress  spatial  variations.  However  a 
number  of  first  order  models  have  been  constructed  and  are  in  use.  They 
appear  to  be  in  rather  surprisingly  good  agreement  with  observations  and 
complex  numerical  models  of  earthquakes. 


2 


The  source  theory  coupled  with  rather  general  wave  propagation  theory 


has  been  used  to  generate  theoretical  and  predictions  for  a wide  class 


of  earthquakes.  These  results,  when  compared  with  comparable  results 


for  explosion  sources,  provide  a detailed  explanation  of  the  mj^ 


discriminant.  These  results  are  discussed  and  illustrated  in  sections  111 


and  IV  of  this  report. 


We  have  also  devised  a variable  frequency  magnitude  discriminate  employing 


magnitude  measurements  at  particular  fequencies,  for  both  body  and  surface 


waves,  with  a pair  of  such  measurements  used  in  the  same  numner  as  the 


(m^,  M^)  pair  for  discrimination  of  explosions  and  earthquakes.  A high 


frequency  body  wave  discriminate,  using  an  event  nuignitude  pair  measured  near 


.4  Hz  and  near  3 Hz  has  been  studied  extensively  using  hundreds  of  actual 


events  from  Eurasia.  We  find  this  discriminate  to  be  superior  to  the 


m^  ^ approach  and  to  be  applicable  over  a wider  event  magnitude  range- 


The  method  is  also  less  costly  to  implement  and  has  been  computer 


automated  so  that  it  is  easy  and  cheap  to  apply.  This  approach  is  discussed 


in  sections  HI  and  IV, 


Several  hundred  events  from  the  Alaska-Aleut Ians  and  Japan-Kamchatka 


regions  were  used  to  infer  the  tectonic  stress  levels  within  the  lithosphere. 


The  theoretical  ^ results  were  used  as  the  basis  for  the  Inference 


of  stress  drop  and  final,  maximum,  failure  zone  dimension.  Moderate  sized 


events  were  used  in  this  study  in  order  to  avoid  uncertainties  associated 


with  strongly  variable  prestress  and  stress  drop.  These  latter  problems 


and  considerations  are  discussed  in  the  following  section  II  while  the 


results  of  the  stress  study  are  summarized  in  section  V. 


Much  of  the  work  reported  has  been  carried  on  in  a joint  research  effort 


with  personnel  at  Systems,  Science  and  Software  and  with  J.B.  Minster  at 


I 
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II.  Theoretical  Investigations;  Seismic  Source  Theory  and  Models. 

Relaxation  theory  models  have  been  used  to  provide  a description  of  the 
radiation  field  from  an  earthquake.  The  theory  is  described  by  Archambeau 
(1968,  1972)  and  Minster  (1973).  The  Appendix  A contains  a very  general 
formulation  of  the  theory  that  incorporates  the  dynamical  boundary  conditions 
explicitly  in  a general  Green's  function  formulation  (Archambeau  and  Minster, 
1977).  Minster  (1973)  and  Minster  and  Archambeau  (1977)  have  incorporated 
the  theory  in  computer  codes  which  allow  detailed  computation  of  the  radiated 
field  from  a number  of  models  and  these  have  been  combined  with  body  and 
surface  wave  propagation  methods  to  allow  the  field  to  be  predicted  at  both 
near  and  teleseismic  distances,  in  both  the  time  and  frequency  domains 
(Cherry  et  al.,  1973;  Cherry  et  al.,  1974;  Savino  et  al . , 1975;  Bache  and 
Archambeau,  1976;  Bach  and  Harkrider,  1976).  The  wave  theory  has  been  adapted 
from  a number  of  sources  and  utilizes  generalized  ray  methods  (e.g.,  Gilbert 
and  Helmberger,  1972;  Wiggins  and  Helmberger,  1973)  as  well  as  classical  ray 
theory  and  full  wave  theory  for  both  body  and  surface  waves  (e.g.,  Bach  and 
Harkrider,  1976;  Harkrider,  1964;  and  Harkrider  and  Archambeau,  1976).  The 
wave  propagation  theory  is  reasonably  straight  forward  and  we  use  either  a 
generalized  ray  theory  or  full  wave  theory  when  required — such  as  in  the  near 
field  distance  range  from  the  source,  at  caustics  at  teleseismic  distance 
ranges  and  in  the  computation  of  the  crustal  transfer  functions  at  the  source 
and  receiver.  For  studies  of  some  classes  of  events,  in  particularly 
complex  structural  regions,  we  can  employ  three  dimensional  ray  tracing 

(e.g.,  Engdahl,  1973;  and  Engdahl  and  Lee,  1975). 

Current  relaxation  theory  models  are  also  being  refined  to  include 

effects  such  as  spatially  variable  prestress  (and  stress  drop)  and  variable 
rupture  velocity.  They  are  also  being  extended  to  account  for  the  dynamic 
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boundary  conditions  connecting  these  variables  (Archambeau  and  Minster, 

1977).  In  addition,  source  models  with  ellipsoidal  rupture  zones  and  with 
variable  failure  material  properties  will  be  included,  with  all  these 
complications  being  incorporated  through  the  use  of  combined  numerical- 
analytical  methods.  These  developments  were  partially  completed  during 
the  contract  period  of  this  report  with  the  results  to  be  incorporated 
as  the  models  are  perfected.  It  is  felt  that  x i view  of  the  clear  evidence 
for  complex  failure  behavior,  such  as  variable  stress  drops  along  the 
failure  "plane"  (e.g.,  Jungels,  1973,  Jungels  and  Frazier,  1973;  Hanks,  1974) 
as  well  as  observations  of  "multiple  events"  which  is  another  (extreme) 
manifestation  of  variable  stress  drop;  then  it  is  important  to  obtain 
"higher  order"  models  so  that  we  can  investigate  and  assess  the  importance 
of  these  effects  in  detail  (see  also  Archambeau,  1975,  for  additional 
comments) . 

In  this  regard,  however,  we  have  compared  predictions  of  the  currently 
used  relaxation  theory  models  to  numerically  computed  radiation  fields  from 
relatively  complex  two  and  three  dimensional  failure  models  incorporating 
plasticity  and  a rather  general  failure  condition  (Cherry,  Bache  and  Masso, 
1976  and  Cherry',  Halda  and  Hamilton,  1976).  The  fields  obtained  from  the 
numerical  calculations  were  compared  to  those  predicted  from  the  analytical 
relaxation  theory  models,  for  which  the  same  prestress,  rupture  velocity 
and  maximum  rupture  dimension  were  used,  and  it  was  found  that  the  two 
results  were  essentially  the  same.  This  implies  that  the  analytical 
models,  even  though  idealized  and  employing  a failure  zone  of  circular 
cross-section,  are  good  approximations  to  failure  situations  in  which  the 
failure  region  is  both  narrow  and  subject  to  nonlinear  flow  laws.  This 


agreement  also  suggests  that  the  radiated  field  is  insensitive  to  the 
thickness  of  the  failure  zone  used  in  the  models,  so  long  as  the  overall 
surface  area  of  the  failure  region  is  comparable  to  that  for  the  (expected) 
thin  zone.  Further  the  results  imply  that  the  radiated  field  is  also 
insensitive  to  any  localized  plastic  flow  occurring  near  the  rupture  front 
prior  to  actual  failure,  so  long  as  the  failure  rate  used  for  the  analytical 
model  is  comparable  to  the  rate  of  actual  failure  zone  growth.  The  numerical 
modeling  did  not,  however.  Include  spatial  variations  of  prestress,  so  no 
detailed  information  on  the  effect  of  pres.'ress  concentration  was  obtained; 
nor  could  we  therefore  make  comparisons  to  the  approximations  used  in  current 
relaxation  source  theory  models  to  account  for  stress  localization  effects 
(i.e.,  to  evaluate  the  "R^  approximation"  as  introduced  by  Archambeau, 

1968.)  However,  additional  work  in  this  area  is  in  progress  and  we  hope  to 
evaluate  its  influence  on  the  radiated  field  in  general  circumstances. 

The  spectral  complications  of  the  observed  seismic  field  from 
earthquakes  arc  Indicated  by  the  spectra  shown  in  Figure  (1).  The  data 
were  recorded  in  the  near  field  distance  range  by  the  University  of  V.'ashLngto 

and  displacement  spectra  for  several  eartliqnakcs  are  shown.  The  detectorr: 
were  long  period  vertical  seismometers,  and  the  noise  at  high  frequency  is 
rather  high. 

It  is  apparent  that  all  the  spectra  are  quite  strongly  peaked  and 
that  at  low  frequency  the  spectra  increase  as  f”*^,  with  n ^ 1.  The  latter 
effect  is  expected  to  be  f , due  to  the  static  offset  in  displacement 
from  "near  field"  terms,  and  the  observed  variation  may  be  a combination  of 
the  expected  behavior  and  noise  (of  various  origins).  The  peak  in  tiie 
spectra  may  be  due  to  interference  due  to  wave  propagation  effects  giving 
rise  to  cancellation  and  a spectral  "hole"  (between  -1.5  and  -1.0  on  the 
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log  frequency  axis  in  the  figure)  or  it  may  be  due  to  spatial  variations  in 
the  initial  tectonic  stress  field  whicli  gives  rise  to  a source  related 
interference  effect  and  a peak  in  the  spectrum.  Tlie  latter  seems  more  likely 
since  it  occurs  for  all  the  events  as  well  as  for  other  receivtfrs  (see 
Figure  2). 

A fit  to  an  observed  radiation  spectrum  from  an  earthquake  is  shown 
in  Figure  2.  The  distance  factor  (which  is  used  in  the  approximate  theory 
to  represent  the  characteristic  spatial  dimension  of  a zone  of  high  stress 
concentration  in  the  vicinity  of  the  failure  volume)  accounts  for  the  "hole" 
in  the  spectrum  near  .1  Hz,  although  its  true  effect  is  to  cause  a peak  in 
the  displacement  spectrum  so  that  the  "hole"  is  actually  due  to  a higher 
frequency  peak  followed  by  the  near  field  caused  increase  in  the  spectral 
level  at  much  lower  frequencies  (below  .1  Hz,  in  this  example).  In  any 
case  the  fault  dimensions  and  rupture  velocity  are  reasonable  for  this 
event  and  a prestrain  level  (e)  of  5 x 10  ^ implies  a stress  drop  of 
about  150  bars  for  a rigidity  change  before  and  after  failure  of  about 
3 X 10^^  dynes/cm^. 

It  should  be  noted  that  while  the  spectral  observations  illustrated 
in  Figures  1 and  2 indicated  complications  at  the  lower  frequencies  giving 
rise  to  a spectral  peak  (or  a spectral  minimum  depending  on  how  one  wishes 
to  describe  the  spectra),  this  effect  is  not  always  observed  for  earthquakes  - 
at  least  over  a fairly  wide  range  of  frequencies  for  some  events.  Tucker 
and  Brune  (1973)  for  example  have  obtained  S wave  spectra  for  a large  number 
of  small  aftershock  events  which  are  nearly  flat  over  almost  a decade  of 
frequency  below  the  corner  frequency.  In  the  present  context  this  would 
imply  that  the  spatial  extent  of  the  zone  of  high  prestress  was  large.  In 
fact,  the  spectral  "peak"  is  predicted  to  be  broad  (and  quite  flat)  and  only 
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Figure  2.  Comparison  of  a theoretically  predicted  cartlKiuake  spectrum  with  an 
observed  earthquake  spectrum  (solid  circles).  The  spectra  are 
composites  of  botli  P and  S waves  from  tiie  event  and  contain  both  near 
and  far  field  spectral  components  (i.c.,  the  total  field  spectrum 
of  the  event).  The  spectrum  was  observed  at  a point  only  14  km  from 
f'  ttie  hypocenter  of  the  eartl\quake.  I'arris  Ranch  earthquake  (California) 

of  27  Oct.  1969.  (Data  from  Johnson  and  McCvilly,  1974.) 
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sharply  peaked  when  the  dimension  (R^)  of  tlie  high  stress  zone  Is  comparable 

to  the  final  fault  dimension.  Tiie  width  of  the  spectral  maximum  is  given 

to  first  order  by  the  difference  in  tlie  corner  frequency  f V /i.  and  the 

frequency  f ~ V'  /R,.  (for  sliear  waves  where  V..  is  tlie  medium  sliear  velocity). 

S S S 

Thus  if  the  dimension  of  the  high  prestress  zone  is  large,  tlien  the  width 
of  the  spectral  peak  can  be  large.  For  small  events  (i.e.,  where  L is  small) 
this  is  particularly  likely  to  be  the  case. 

Thus,  spectral  complications  (e.g.,  spectral  maxima  and  minima) 
at  frequencies  below  the  corner  frequency  may  occur,  but  they  may  be 
less  likely  to  occur  in  the  observational  band  width  used  for  small  earthquake 
detection  if  they  are  associated  with  stress  concentration  effects.  If  they 
are  associated  with  structural  induced  wave  interference  effects  - and  at 
higher  frequencies  they  definitely  are  - then  we  can  avoid  any  biasing  effects 
in  spectral  magnitude  measutements  by  proper  averaging  over  several  observations 
at  different  distances  and  aximuths  and  by  use  of  several  spectral  estimates 
and  smoothing. 

However,  there  clearly  will  be  events  observed  with  spectra  like 
those  of  Figures  1 and  2 and  such  simple  procedures  as  averaging  and  smoothing 
may  not  remove  the  bias  effect  of  spectral  magnitude  measurements  made  in 
the  spectral  minimum,  when  the  minimum  is  a result  of  a spatial  stress 
variation  effect.  Further,  this  minimum  may  occur  for  both  P and  S waves 
j and  will  also  then  affect  ( reduce)  the  ordinary  surface  wave  magnitude 

value,  M^,  for  such  events.  We  therefore  seek  to  avoid  making  any  measurement 
in  this  zone  for  the  class  of  small  earthquakes  we  Intend  to  use  for  stress 
estimation,  whether  or  not  the  effect  occurs.  We  can  also  use  spectral 
magnitudes  for  discrimination  that  are  unlikely  to  be  affected.  This  is 


11 


not  difficult  to  do  if  the  characteristic  dimension  of  the  uniformlv  high 
initial  stress  zone  is  of  the  order  of  10  to  20  km  and  larger  in  the  vicinity 
of  earthquakes.  In  particular.  Figure  3 shows  tlieoretical  spectra  for 
earthquakes  with  small  nuiximum  dimensions  (from  1 km  up  to  10  km)  where  the 
characteristic  dimension  for  the  stress  concentration  varies  with  the  fault 
dimension  from  10  km  to  100  km.  From  this  figure  we  see  that  for  small 
earthquakes,  if  we  use  a low  frequency  magnitude  estimate  near  .3  Hz  and  a 
high  frequency  magnitude  estimate  near  5 Hz,  we  can  be  quite  confident  that 
the  low  frequency  estimate  will  lie  on  or  near  the  broad  nuaximum  of  the  P 
wave  spectrum  while  the  high  frequency  estimate  will  be  above  the  corner 
frequency  on  the  declining  high  frequency  slope.  In  this  case,  the  only 
parameters  controlling  the  spectrum  in  this  range  are  the  rupture  dimensions, 
the  stress  drop  and  the  rupture  velocity  and  in  terms  of  the  measured  spectral 
magnitudes  used  for  discrimination  and  source  parameter  estimates  the  events 
appear  as  normal  earthquakes.  Normally  for  stress  estimates  the  rupture 
velocity  is  assumed  to  be  about  .8  of  the  local  shear  velocity,  so  the  remaining 
two  parameters  affecting  the  measured  spectral  magnitude  are  the  rupture 
length  and  stress  drop.  These  event  parameters  can  then  be  determined  using 
pairs  of  spectral  magnitude  estimates  - as  will  be  discussed  in  more  detail 
in  the  following  sections.  Further  use  of  the  "high"  frequency  spectral 
magnitudes  for  discrimination  avoids  or  minimizes  the  possible  problem  of  a 
spectral  minimum  at  low  frequency  causing  a depression  of  the  low  frequency 
magnitude  estimate  (such  as  could  occur  for  an  measurement).  Such 
depression  could,  when  low  and  high  frequency  magnitudes  are  used,  cause 
merging  of  the  earthquake  population  and  the  explosion  population.  Hence 
use  of  a high  frequency  body  wave  spectral  magnitude  discriminant  appears 
to  be  better  approach  to  explosion-earthquake  discrimination  then  the 
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Figure  3,  Scaling  of  the  epectrun  for  different  fault  lengths, 
v;ith  the  ratio  of  the  semi-major  axis  length  to 
major  axis  length  of  the  fault  held  fixed. 
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method  for  this  reason,  as  well  as  because  It  is  simpler  and  has  a 
greater  magnitude  range  of  applicability. 

If  the  spatial  stress  concentration  has  a characteristic  dimension  of  the 
order  of  or  greater  than  50  to  60  km  for  all  events  tlien  tmignitudes  would 
not  be  affected  by  any  spectral  minimum  associated  with  the  prestress. 

It  may  be  the  case  that  most  events  satisfy  this  condition,  however  it  is 
quite  possible  tliat  some  small  events  are  the  result  of  failure  within 
small  zones  of  high  stress  concentration.  In  general  then, it  is  safest 
to  use  as  high  a frequency  as  possible  for  the  low  frequency  magnitude 

used  in  discrimination  with  the  variable  frequency  magnitude  . 

i 

approach. 

The  remaining  theoretical  problems  relating  to  discrimination  and  ' 

source  parameter  estimates  are  largely  connected  with  the  stress  concentration  ' 

effects  as  has  already  been  observed.  There  is  still  considerable  controversy 

I 

(see  for  example  Archambeau,  1975)  about  the  manifestations  of  variable  ! 

I 

stress  drop  and  prestress  in  the  radiated  seismic  field  from  an  earthquake. 

It  has  therefore  been  our  approach  in  this  study  to  design  methods  of  stress  | 

estimation  and  event  discrimination  that  avoid  use  of  data  that  may  be  j 

effected  by  the  stress  spatial  variability,  while  simultaneously  attempting  i 

to  resolve  the  theoretical  questions  in  a rigorous  fashion.  The  later  effort  , 

seems  to  be  quite  close  to  a final  solution,  but  will  require  somewhat  more 

effort.  We  think,  at  the  moment  at  least,  that  the  approximate  method  involving 

the  use  of  a characteristic  dimension  (R  ) for  a stress  concentration 

s 

(E.g.,  Archambeau,  1968)  is  roughly  correct  in  terms  of  the  predicted 
effects  on  the  radiation  field,  but  of  course,  this  cannot  be  rigorously 

I 

supported  now.  I 
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In  applications  of  the  theory  to  predictions  of  and  for  use  in 
discrimination  and  stress  estimates,  we  have  assumed  that  the  prestress  zone 
is  uniform  over  a region  of  50-60  km,  so  that  the  value  is  not  affected 
by  any  stress  concentration  effects  and  also  have  confined  ourselves  to  the 
analysis  of  moderate  sized  earthquakes  where  this  is  very  like  to  be  true. 

A major  part  of  our  theoretical  work  has  been  directd  to  the  formulation 
of  relaxation  source  theory  in  a rigorous  and  general  manner  that  explicitly 
accounts  for  the  dynamical  boundary  conditions  on  an  expanding  failure 
boundary.  We  have  also  introduced  the  idea  of  a generalized  phase  change 
and  an  associated  internal  energy  change  in  material  parameters  as  a means 
of  quantitatively  describing  the  energetics  of  a failure  process.  This 
work  is  described  in  detail  in  the  Appendix  A of  this  report. 


15 


r 

III.  Sumnmry  of  ^nd  ">^(0  Mflhods  of  Stress  Estimation  and  Kvent 

Discrimination. 

A basis  for  discrimination  and  identification  of  explosions  and  earthquakes, 
and  estimation  of  source  properties,  in  particular  stress  drop  of  earthquakes 
and  energy  yield  of  explosions,  can  be  established  througli  the  use  of: 

(1)  A realistic  and  accurate  theoretical  repre.sentat ion  of  the  source 
radiation  field. 

(2)  A reasonably  accurate  and  detailed  knowledge  of  the  crust  and 
upper  mantle  elastic  and  anelastic  properties. 

(3)  The  capability  to  generate  theoretical  spectra,  and  seismic  time 
series  for  at  least  the  major  wave  types  comprising  the  seismic  field  observed. 

(4)  A large  data  base  containing  multiple  seismic  observations  of 
small  to  moderate  sized  events. 

(5)  An  efficient,  fast  and  accurate  process  for  estimating  the  important 
observational  parameters  to  be  used  for  event  discrimination  and  in  the 
estimation  of  source  parameters  - in  particular  stress  for  a large  event  set. 

The  first  three  of  these  elements  are  used  to  establish  a theoretical  basis 
for  event  discrimination  and  source  parameter  inversion  and  the  final  two  are 
required  to  verify  the  theoretical  approach  and  to  implement  it.  In  this 
section  we  will  discuss  the  use  of  our  theoretical  modeling  ability  in 
establishing  methods  of  event  discrimination  and  event  parameter  determinations, 
in  particular  the  stress  level  changes  associated  with  earthquakes.  In 
the  previous  section  we  have  discussed  the  theoretical  source  models 
that  will  be  used.  In  the  following  sections  we  will  describe  the 
signal  analysis  methods  used  in  generating  appropriate  observational  data 
and  then  illustrate  the  approach  through  applications  giving  source  parameter 


estimates  (section  V)  and  event  discrimination  (section  VI).  The  use  of 
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the  observed  data  is  also,  of  course,  designed  to  test  the  accuracy  and 
reliability  of  the  theoretical  methods  devised  for  discrimination  and  source 
parameter  estimation.  (0  ir  verification  criteria  is  based  on  the  internal 
consistency  of  the  source  parameter  estimates  as  well  as  upon  the  predictive 
reliability  of  the  results.) 

Theoretical  m^^  and  M_  Magnitudes  for  Earthquakes, 

The  source  models  described  and  the  ability  to  predict  the  teleseismic 
radiation  from  such  models  in  the  earth  have  been  used  to  generate 
theoretical  or  "synthetic"  seismograms  for  earthquakes  of  various  types. 

In  particular,  theoretical  seismograms  have  been  generated  for  earth- 
quakes at  variable  depths,  in  the  range  0-50  km,  with  variable  length, 
orientation  (l.e. , strike-slip,  dip-slip,  normal  and  thrust  types) 
and  with  variable  stress  drops.  The  predicted  radiation  fields — consisting 
of  compressional  body  waves  and  surface  waves — from  this  variety  of 
earthquakes  have  been  computed  over  a range  of  teleseismic  distances 
and  at  numerous  azimuths.  These  detailed  predictions  have  then  been 
used  to  obtain  m^^  and  values — measured  from  the  synthetic  data 

and  averaged  by  the  identical  procedures  used  in  routine  observational  practice. 
The  details  of  this  computational  program  are  described  by  Archambeau 
et  al.  1974,  Archambeau  (1974,  1975b),  Bache  and  Archambeau  (1976) 

Savino  and  Archambeau  (1976)  and  Archambeau  (1976), 

Figures  ( 4 ) through  ( 6)  are  examples  of  synthetic  seismograms 
for  a dip  slip  earthquake  with  hypocentral  depth  of  10  km  and  stress 
drop  of  500  bars.  The  rupture  velocity  for  the  event  was  taken  to  be 
.9v  , where  v is  the  shear  velocity  at  the  10  km  depth  in  the  earth 

S 8 

model  used  for  the  structure  at  the  source  location.  In  this  example, 
a continental  tectonic  model  was  used  (CIT  109  , Archambeau  et  al. , 1969)  . 
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Figure  ( A ) illustrates  some  of  the  complications  that  arise  from  the 

rapid  velocity  transition  zones  in  the  upper  mantle  and  fronr  the 

rather  complicated  character  of  the  signals  generated  within  the  source 

region.  Both  direct  compressional  waves  and  converted  waves  (eg.  shear 

to  compressional  reflections)  are  represented  in  the  wave  trains  shown. 

Figure  ( 5 ) shows  predicted  Rayleigh  waves  from  several  earthquakes  of 

the  type  shown  in  Figure  ( A ) , but  with  failure  zone  length  ranging  from 

1 to  20  km.  Figure  (6^  shows  the  Love  waves  from  this  same  set  of  events. 

The  measurements  of  and  are  of  course  made  from  parts  of  the 

compressional  and  Rayleigh  wave  trains  illustrated.  Generally,  the  body 

wave  measurement  is  made  from  the  largest  cycle  in  the  first  three  cycles 

of  the  compressional  wave  train.  For  shallow  earthquakes  this  often 
results  in  a measurement  of  the  surface  reflected  shear  to  compressional 

wave  rather  than  a measurement  of  the  direct  compressional  wave 

from  the  source.  Since  this  may  occur  in  observational  practice  it  is 

reflected  in  our  procedure  of  m^  measurement  as  well. 

In  this  regard.  Figure  (7)  illustrates  the  variation  in  predicted 

teleseismic  recordings  for  earthquakes  with  variable  hypocentral  depths, 

indicated  by  H,  ranging  from  10  to  A5  km,  in  the  figure.  The  "phase" 

used  to  compute  the  m^  value  associated  with  the  individual  recordings 

is  indicated  by  a bar  or  slash  at  the  zero  crossing  between  the  peak 

to  peak  values  used  in  the  computation.  Clearly  the  large  secondary  phase 

is  a converted  signal — specifically  the  free  surface  converted  shear  to 

compressional  wave — and  for  the  event  at  H = 10  km,  this  phase  would 

be  used  to  compute  the  m value  using  the  standard  conventions  applied 

b 

to  magnitude  computations.  This  phase  might  also  be  used  for  the  15  km 


depth  event  as  well. 
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Figure  7:  Theoretical  seismograms  at  an  epicentral  distance  of  .050 


of  this  phase. 
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Aside  from  illustrating  the  variation  in  the  seismic  observations 
as  functions  of  source  depth.  Figure  (/)  provides  an  indication  of 
the  variations  to  be  expected  as  a function  of  azimuth  from  tlie  event. 
Since  the  two  samples  in  azlnutli  are  separated  by  90°,  we  observe  that 
one  time  series  is  close  to  being  tlie  negative  of  the  otiier;  that  is, 
to  first  order  the  corresponding  spectra  are  uniformly  different  in 
phase  by  180°  These  reversals  in  sign,  for  the  "first  motion" 
section  of  the  signal  train,  is  of  course  the  basis  for  the  commonly 
used  fault  plane  solutions.  In  addition,  however,  there  are  differences 
in  wave  form  and  frequency  content  in  the  synthetics  at  the  different 
azimuths,  which  result  in  variations  in  the  effective  period  T 
measured  for  the  so  called  "b  phase"  and  its  amplitude  and  consequently 
in  the  values  computed  at  the  two  azimuths.  These  effects  are 

due,  principally,  to  rupture  propagation  effects. 

The  variation  in  the  theoretical  seismograms  with  azimuth  are 
also  shown  in  Figure  (8)  for  two  different  "types"  of  earthquake, 
that  is  for  earthquakes  oriented  in  different  directions  with  rupture 
propagation  and  shear  stress  fields  correspondingly  oriented  in 
different  directions.  Here  the  failure  "plane"  for  the "strike-slip" 
event  is  perpendicular  to  the  earth's  free  surface  while  the  sense  of 
rupture  growth  is  parallel  to  the  surface  with  an  initial  shear  stress 
field  oriented  such  that  the  shear  offset  along  the  failure  zone  is 
parallel  to  the  free  surface.  The  failure  plane  for  the  "dip-slip" 
event  is  at  45°  with  respect  to  the  free  surface  and  with  rupture 
propagation  down  the  plane  away  f.om  the  free  surface  , with  an  offset 

due  to  the  prestress  orientation  which  is  at  45°  to  the  surface. 
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Figure  8:  Theoretical  seismograms  for  several  azimuths  for  an  earthquake  source  at 
two  orientations,  strike-slip  and  dip-slip.  The  source  depth  is  25  km  and  the 
epicentral  distance  is  4000  km. 
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The  theoretical  compress ional  wave  time  series  for  these  two  events 
are  quite  different  as  would  be  expected.  Both  events  have  the  same 
rupture  length  (1  km) , the  same  rupture  velocity  magnitude  and  the  same 
stress  drop  (100  bars).  Hence,  the  differences  In  the  amplitudes  of 
the  b phase  used  to  compute  m^,  for  example,  for  the  two  source 
types  is  due  entirely  to  the  event  orientation.  Further,  the  strike- 
slip  event  shows  strong  azimuth  variation,  with  change  in  sign  of  the 
"first  motion"  of  the  signal  train  occurring  between  75“  and  105“, 
while  the  dip-slip  event  shows  little  variation  in  wave  form  with 
azimuth  but  substantial  variation  in  amplitude.  These  effects  are  due 
to  the  radiation  pattern  from  the  source  in  combination  with  interference 
ef..ects  due  to  multiple  reflections  within  the  medium  near  the  source 
(principally)  and  at  the  receiver. 

VJhile  the  amplitude  and  wave  form  variations  with  azimuth  and  event 
type  are  seen  to  be  rather  large  effects,  the  variation  in  the  average 
amplitude  of  the  P wave  train  as  a function  of  depth  is,  as  might  be 
expected,  not  large.  This  is  illustrated  by  Figure  (9).  Also  shown 
is  the  change  in  character  of  the  wave  train  for  two  distances,  one  at 
2500  km,  where  upper  mantle  velocity  gradients  produce  multiple  direct 
arrivals  and  the  other,  at  4000  km,  where  these  complications  are  absent. 
The  seismograms  at  2500  km  show  quite  clearly  that  the  later  arriving 
energy,  corresponding  to  direct  compressional  waves  resulting  from  the 
sharp  velocity  gradient  near  450  km  in  the  mantle,  are  considerably 
larger  than  the  first  arriving  compressional  wave  signal.  In  some 
cases  at  least,  this  signal  would  be  used  to  compute  the  in  magnitude. 
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Figure  9:  Theoretical  seismograms  for  five  depths  and  two  epicentral  distances  for  the 
relaxation  source  model  of  Archambeau  (1968) < The  faulting  is  strike-slip  and  the 
azimuth  is  30°  from  the  strike. 
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If  we  use  the  current  definition  for  m and  assume  that  it  is  followed 

b 

by  those  observers  that  report  la^  values,  then  since  the  synthetic 
seismograms  contain  all  the  essential  complications  of  the  observations, 
we  can  use  the  theoretical  time  series  to  obtain  m^  values  accounting 
for  these  complications.  Employing  the  observational  definition  of  m^ 
throughout,  individual  station  values  at  the  appropriate  distances  and 
azimuths  relative  to  the  source  may  then  be  averaged  to  give  a theoretical 
m^  for  the  event.  The  same  procedure  can  be  used  to  generate  theoretical 

M values, 
s 

Using  theoretical  results  such  as  those  shoxm  in  Figures  ( 4 ) through 
(9),  we  therefore  can  generate  the  required  m^  - results  for  a 
variety  of  seismic  events  at  various  depths.  Figure  (10)  shows  m^  - 
results  for  dip-slip  earthquakes,  with  fault  length  and  stress  drop  variable 
while  the  rupture  velocity  and  hypocentral  depths  are  held  constant,  for 
all  events.  Thus  the  lines  labeled  .01  k bar  to  1 k bar  are  the  - M 
loci  of  events  of  variable  dimension  and  fixed  stress  drop,  while  the 
second  set  of  lines  define  loci  of  events  with  fixed  dimension  but  with 
variable  stress  drop.  These  curves  obviously  define  a coordinate  grid 
in  the  m^  - plane  which  can  be  used  to  interpret  m^^  and  M ^ values 
in  terms  of  rupture  dimension  and  stress  drop.  Such  theoretical  results 
can  of  course  be  generated  for  different  event  types,  at  different 
depths  and  with  a different  (fixed)  value  of  rupture  velocity  if  necessary. 
This  set  of  theoretical  results  then  provides  the  framework  for  the 
interpretation  of  m^  - observations  in  terms  of  stress  drop  and 
failure  zone  dimension.  Further,  if  we  assume  that  the  largest  stress 
drops  observed  for  a region  are  "total"  stress  drops,  corresponding  to 
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Figure  10:  Theorot ieal  roj^-Mg  curves  for  dip-slip  earthquakes  at  10  km  depths  all  with 
the  same  rupture  velocity  \^=  but  having  variable  maximum  dimension  and  stress  drop 
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a complete  loss  of  shear  strength  within  the  failure  zone  during  failure, 
then  the  largest  stress  drop  observations  for  a given  local  region  define  the 
initial  stress  field. 

Figure  (11)  illustrates  the  use  of  the  theoretical  curves  for  observed 
shallow,  dip-slip  events  in  the  Alaskan-Aleutian  region.  (The  lines 
of  constant  failure  dimension  have  been  omitted  for  clarity,  but  are 
parallel  to  the  shown.)  For  these  dip-slip  events,  the 

average  stress  drop  is  near  100  bars. 

The  dotted  line  labeled  "NTS"  is  an  - M line"  for  explosions, 
with  explosion  yield  variable  along  this  line.  (A  set  of  such  explo- 
sion lines  exist,  corresponding  to  explosions  in  different  materials.) 

Theoretical  nij^  - relations,  such  as  are  shown  in  Figures  (10)  and 
(11),  are  therefore  a basis  for  stress  level  estimates  as  well  as  for  the 
discrimination  of  earthquakes  and  explosions. 

Theoretical  Variable  Frequqncy  Magnitudes. 

As  either  an  alternate  or  a supplement  to  the  "m^  - method"  for 

discrimination  and  estimating  tectonic  stress,  we  will  use  what  can  be  termed  a variable 

frequency  magnitude  method,  or  method".  This  simply  refers  to 

the  use  of  body  wave  magnitudes  measured  from  output  of  ultra  narrow  band 

filters  at  two  (or  more)  frequencies.  Otherwise,  the  inference  of 

source  parameters  is  essentially  the  same  as  described  for  the  m^  - M 

method,  but  with  m^^  measured  at  two  frequencies  used  in  place  of  the 

(nijjjMs)  pair.  Here  we  use  m^(f)  measured  at  a moderately  low 

frequency — such  as  .1  or  .2  cps — in  place  of  M . The  role  of  the 

s 

traditional  m^  measurement  is  replaced  by  in^(f)  measured  at  a high 
frequency — in  the  range  1.0  to  3.0  cps.  usually. 
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Figure  11:  Theoretical  and  observed  data  for  Alaskan-Aleutlan  Arc  region 

earthquakes  in  the  depth  range  from  0 to  fo  km.  The  curves  shown  on  the  figure  are 
parameterized  with  respect  to  the  stress  drop  with  the  rupture  dimension  variable  along 
each  of  the  curves. 
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This  approach  was  designed  by  Archambeau  et  al.  (1974)  as  a means 
of  discriminating  shallow  earthquakes  from  explosions  and  has  been 
applied  extensively  for  this  purpose.  Savino  and  Archambeau  (1976), 


Archambeau  (1975b),  Bache  et  al.  (1975)  and  Cherry  et  al.  (1974b), 
show  numerous  examples  of  the  generation  and  use  of  in,  (f)  data. 


Figure  (12)  shows  an  example  of  the  narrow  band  filtering  applied 
to  an  observed  explosion  and  an  earthquake.  We  observe  that  peaks  in 
the  narrow  band  filter  output  can  be  correlated  in  time  with  the  onset 
of  the  original  compressional  wave  signal.  The  amplitudes  of  these  time 


correlated  peaks  in  the  filter  output  may  be  used  to  compute  the  m (f) 

b 


values,  and  give  magnitudes  to  be  associated  with  the  frequency  defined 
by  the  center  frequency  of  the  filter.  Use  of  narrow  band  filter  output 
amplitude  data  (as  compared  to  spectral  amplitudes  from  the  complete 
Fourier  spectrum  analysis  of  the  entire  compressional  wave  train)  is  there- 
fore advantageous  in  that  time  information  is  retained  in  the  filter  output, 
so  that  the  amplitude  measured  can  be  related  to  a particular  transient 
signal  arrival.  (This  is  not  the  case  for  a Fourier  analysis  of  a large 
time  segment  containing  several  transient  signals.)  \7hile  the  time  infor- 
mation is  retained,  the  amplitude  measurement  can  still  be  related  to  a parti 
ular  frequency  component  of  the  signal.  Here,  of  course,  there  is  a "trade- 
off" between  time  and  frequency  resolution,  which  may  be  stated  as  AfAt  > 1, 
where  At  is  the  uncertainty  in  the  time  of  the  frequency  arrival  (arrival  tin 
of  the  energy  at  the  frequency  f)  and  with  Af  the  uncertainty  in  the 
frequency  determination.  Nevertheless  we  can  allow  Af  for  each  amplitude 
estimate  to  be  finite  (it  is  proportional  to  the  filter  width  at  the 
half  power  points)  and  still  obtain  a very  good  estimate  of  spectral 


Figure  12:  Application  of  a "comb"  of  narrow  band  filters  to  seismic  events  to  obtain 
spectral  amplitudes  and  time  of  energy  arrival  for  several  frequencies  In  the  seismic- 
band.  An  explosion  and  deep  earthquake  are  Illustrated.  The  amplitudes  obtained  are 
used  to  compute  variable  frequency  magnitudes  mj^Cf)  for  P waves. 
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changes  in  the  signal  with  frequency.  Then  with  finite  Af,  we  see 

that  At  will  also  be  finite  so  that  we  can  retain  time  resolution 

and  thereby  be  able  to  associate  the  measured  amplitude  with  a particular 

part  of  the  rather  complex  series  of  signal  pulses  comprising  the  wave 
train,  as  was  the  case  with  the  classical  m^^  measurement.  By  this 

process,  we  "buy"  frequency  resolution  in  the  determination  of 

compressional  wave  magnitudes  which  we  may  associate  with  the  first 

arriving  P wave  transient.  In  essence  then,  this  is  a simple  method  of 

obtaining  spectral  amplitude  (or  magnitude)  estimates  of  the  direct 

compressional  wave  signal  from  seismic  events  and  a method  which  is 

easily  automated,  so  that  we  may  deal  with  large  amounts  of  data. 

Consequently,  we  may  obtain  an  estimate  of  the  source  amplitude  spectrum 

at  high  and  low  frequencies  from  the  P wave  alone , without  use  of 

the  surface  wave  magnitude  measurement  as  a low  frequency  measure,  and 

therefore  can  obtain  the  essential  information  required  for  discrimination  and 
stress  estimates  much  more  simply  than  with  the  m^^  - method.  Further, 
this  procedure  can  be  applied  to  very  small  events  and  at  any  distance 
range.  Finally,  it  is  a more  precisely  defined  measurement  than  the 
usual  m^  and  measurements,  so  that  it  offers  the  opportunity  of  greater 

precision  in  the  source  parameter  determinations  and  wider  separation  of 
explosion-earthquake  populations  for  use  in  discrimination  work. 

Narrow  band  filtering  of  the  theoretical  time  series  can  therefore  be  used 
to  define  a theoretical  basis  for  discrimination  and  the  estimation  of  seismic 
event  parameters.  Clearly  the  filtering  operations  can  be  applied  to 
either  the  body  or  surface  wave  signals  to  generate  spectral  magnitudes 
mj^(f)  and  M^(f).  From  these  results,  a variety  of  theoretical  M (f) 
versus  relationships  can  be  generated.  These  relations  are  defined  in 
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a space  in  which  seismic  events  are  coordinate  points  determined  by  their 
ra^^Cf)  and  values,  or  values  at  two  (widely)  separated  frequencies. 

Figure  (l3)  shox^s  a number  of  possibilities  for  one  type  of  event 

at  fixed  depths  in  the  earth.  In  this  figure,  spectral  surface  wave 

*“  R L 

magnitudes  at  .05  cps  for  both  Rayleigh  (M  ) and  Love  waves  (M  ) 

s s 

have  been  plotted  against  1*5  cps  and  2.5  cps,  (labeled 

(1.5)  and  (2.5)  in  the  plots).  Theoretically  predicted  loci 
for  earthquakes  of  variable  fault  dimension,  fixed  stress  drop  (500 
bars)  and  fixed  rupture  velocity  are  shown.  In  the  first  figure  at 
the  upper  left,  loci  for  events  with  100  bar  stress  drops  are  also  shown 
and  indicates  the  typical  change  in  position  of  the  event  loci  in  these 
planes  when  the  stress  drop  changes.  The  general  character  of  the 
loci  are  quite  similar  to  the  "nij^  - lines"  previously  described  and 
and  illustrated  in  Figures  (10)  and  (11), 

In  addition  to  theoretical  earthquake  lines  shown  in  Figure  (13), 
theoretical  points  for  explosions  are  also  shown.  They  appear  in  such 
a plot  as  extremely  high  stress  drop  earthquakes  , when  viewed  in  terms 
of  the  theoretical  earthquake  loci.  The  situation  is  again  similar  to 
that  for  the  m^  - relations.  It  is  obvious  from  these  plots  how 
these  magnitude  measurements  may  be  used  to  discriminate  explosions  and 
earthquakes,  since  the  two  types  of  events  fall  into  separate  populations 
in  these  "spaces". 

The  lower  set  of  four  plots  are  of  greatest  interest  for 
applications  and  shown  the  loci  of  events  of  variable  dimension  and 
fixed  500  bar  stress  drop  in  the  plane.  Essentially  parallel 

curves  in  this  plane  result  from  these  same  sized  events  with  different 
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stress  drops.  Those  with  higher  stress  drops  lying  to  the  right,  closer 
to  the  explosion  population,  while  those  with  lower  stress  drops  define 
curves  to  the  left.  These  stress  variable  "coordinate  lines"  are 

quite  widely  spaced  and  allow  stress  drop  estimates  to  be  obtained  for 
all  events  down  to  those  of  very  small  size. 

Essentially  an  unlimited  number  of  such  plots  could  be  generated, 
since  all  possible  frequencies  can  be  used.  However,  they  are  to  a 
large  degree  redundant  in  that  only  a few  of  them  are  needed  to  obtain 
the  basic  source  property  information  sought — namely  the  stress  level, 
the  failure  zone  dimension  and  rupture  velocity.  If  we  assume  that  the 
rupture  velocity  is  essentially  the  same  for  all  the  events  to  be 
observed,  then  this  parameter  is  effectively  eliminated  from  con- 
sideration and  we  need  only  one  pair  of  variable  frequency  magnitudes 
to  obtain  the  desired  source  properties.  The  magnitudes  required  are 
the  largest  low  frequency  value,  selected  from  several  spanning 

the  low  frequency  spectral  range,  and  the  spectral  magnitude  at  the 
highest  frequency  that  is  well  above  the  noise  level.  Because  the 
spectrum  is  expected  to  be  highly  variable  at  high  frequency,  with 
numerous  "peaks  and  holes",  the  high  frequency  magnitude  value  is  to 
be  obtained  by  using  a straight  line  fit  over  several  frequency 
separated  estimates,  so  that  a smoothed  spectral  value  is  obtained. 

In  any  case,  these  two  magnitudes  would  normally  correspond  to 
spectral  amplitude  estimates  on  either  side  of  the  "corner  frequency" 
the  low  frequency  magnitude  being  obtained  at  the  peak  or  flat  spectral 
level  (whichever  occurs)  and  the  high  frequency  estimate  at  a point 
on  the  rapidly  decreasing  high  frequency  spectral  slope  of  the  radiated 
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compressional  wave  spectrum  from  the  source.  For  a fixed  rupture  velocity 

for  all  events,  the  spectral  shape  is  invariant  and  then  only  these 

two  values  are  needed  to  determine  the  stress  and  failure  zone  dimension. 

In  summary,  these  considerations  show  that  the  spectral  magnitude 
estimates  can  be  used  with  great  flexibility  in  discrimination  studies  and  to 
define  source  parameters.  Further,  in  view  of  inherent  redundancy,  multiple 
estimates  can  provide  several  determinations  of  source  parameters;  with  the 
variations  in  these  estimates  providing  a check  on  the  accuracy  of  the 
determinations.  Similarly,  more  than  a single  discrimination  criteria  can  be 
employed  for  identification  of  explosions  from  a set  of  events.  Clearly  this 
data  can  also  be  used  to  check  the  m^  - based  soruce  parameter  estimates 
and  discrimination  as  well  as  to  extend  the  approach  to  the  very  numerous, 
small  events.  The  latter  capability  is  especially  Important  for  discrimination. 

Finally  and  perhaps  most  Important  of  all  for  the  applications,  is 
the  fact  that  the  data  is  not  strongly  affected  by  source  depth  and 

source  orientation  since,  in  addition  to  weak  or  only  moderate  influence  on 
individual  spectral  magnitudes,  the  magnitudes  pairs  at  the  two  frequencies 
are  both  affected  in  the  same  way  by  source  orientation  and  depth  so  that  the 
effects  essentially  cancel  out  when  the  data  is  viewed  in  the  "mj^(f)  plane" 
defined  by  a pair  of  spectral  magnitudes.  Thus,  in  addition  to  being  only 
strongly  influenced  by  the  prestress  and  rupture  dimension  parameters  so  that 
more  reliable  and  accurate  estimates  of  these  parameters  are  possible, 
this  also  means  that  it  is  not  necessary  to  average  over  observations 
at  different  distances  and  azimuths  and  that  single  station  observations 
can  be  used  in  either  the  near  or  far  field  distance  ranges  to  effect 
discrimination  of  events  and  obtain  source  parameters  estimates.  It  is  therefore 
clear  that  the  "mj^(f)  method"  is  by  far  the  more  flexible  and  accurate 
approach  to  general  source  studies.  Including  discrimination. 
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IV . Signal  Analysis  for  Discrimination  and  Stress  Estimation . 

This  section  is  devoted  to  a description  of  seismic  signal  analysis 
methods  and  the  operational  signal  analysis  computer  based  on  these  methods. 
This  program,  specifically  developed  for  seismic  discrimination  studies,  is 
also  directly  applied  in  source  parameter  studies  - in  particular  to  stress 
estimation. 

Depending  upon  the  availability  of  one  or  more  components  of  recorded 
ground  motion  (vertical  and  two  orthogonal  horizontal  components)  from  seismic 
events,  there  are  various  types  of  filtering  available  in  the  signal  analysis 
program.  Since  the  majority  of  seismic  stations  operate  with  only  a vertical 
component  seismometer,  the  narrow  band  filtering  operation  to  be  described 
first  will  be  most  applicable  for  discrimination  applications  and  for 
determination  of  earthquake  source  parameters  based  on  spectral  observations. 
Whenever  three  component  data  are  available,  both  narrow  band  and  polarization 
filtering  techniques  can  be  applied  in  order  to  isolate  body  phases  for 
spectral  computations  and  event  location  computations. 

The  central  feature  of  the  signal  analysis  program  is  the  use  of  narrow 
band  frequency  filters  to  break  up  or  decompose  a time  series  consisting  of 
signal  plus  noise  into  a set  of  quasi-harmonic  modulated  "signals".  This 
set  of  filtered  signals,  one  for  each  filter  center  frequency,  can  then  be 
used  to  determine  the  energy  arrival  time  (the  group  time,  t ) and  amplitude 

O 

of  the  signal  for  each  center  frequency  by  analysis  of  the  time  modulation 
of  the  filter  outputs.  Further,  both  the  instantaneous  phase  and  frequency 
of  the  individual  filter  outputs,  that  is  the  apparent  phase  and  frequency 
of  the  quasi-harmonic  filter  output  as  a function  of  time,  can  be  determined 
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quite  simply.  Thus,  if  two  components  of  the  same  vector  wave  field  are 
analyzed  in  this  fasion,  the  difference  in  the  instantaneous  phase  between 
the  two  components  is  a direct  measure  of  the  polarization  of  the  wave  field. 

Thus,  the  decomposition  of  the  signal  wave  train,  possibly  composed 
of  many  individual  signal  pulses,  into  quasi-harmonic  signals  provides  the 
means  of  determining  arrival  time,  amplitude  and  wave  polarization,  all  as 
functions  of  frequency.  This  then  is  the  basic  signal  information  that  can 
be  used  to  detect  a given  type  of  signal  in  terms  of  its  dispersion  and 
polarization  characteristics  and  to  obtain  its  spectrum  as  well  as  its  time 
and  amplitude  relation  with  respect  to  other  signals  present  in  a complex 
wave  train. 

Figures  14  and  15  are  two  different  versions  of  the  flow  of  operations 
in  the  existing  computer  program.  Figure  14  provides  a verbal  outline, 
while  Figure  15  summarizes  the  key  mathematical  operations  performed  in 
this  program.  More  detailed  developments  of  the  theory  and  operation  are 
given  in  Bache,  ^£.1.  [1975b]  and  Savino,  ^ al.  [1975]. 

Seismic  data  are  read  into  the  program  in  the  form  of  time  series 
(Figures  14  and  15  ) generally  of  about  500-2000  points  in  length.  The 
data  are  then  optionally  detrended,  demeaned  and  tapered.  The  program 
then  selects  the  smallest  power  of  two  which  is  greater  than  the  number  of 
points  input  and  performs  a discrete  Fourier  transform  using  the  Cooley- 
Tukey  [1965]  algorithm.  Both  the  original  time  series  and  the  spectrum  are 
then  plotted. 

Referring  to  Figures  14  and  15  , the  signal  is  next  filtered  in  the 


frequency  domain  by  multiplication  with  a narrow  band  cusp-shaped  filter. 

This  particular  filter  form  is  selected  to  satisfy  two  goals:  (1)  minimum 
width  in  the  frequency  domain,  and  (2)  maximum  ripple  suppression  in  the  time 
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Flowchart  indicating  principal  mathematical 
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domain.  As  a consequence  of  the  sampling  theorem,  or  uncertainty  principle, 
one  cannot  simultaneously  satisfy  these  two  goals  to  arbitrary  precision. 

The  filter  employed  was  selected  for  its  optimal  time  and  frequency  domain 
characteristics  within  this  basic  limitation. 

Once  the  signal  has  been  narrow  band  filtered,  it  is  corrected  for 
the  appropriate  instrument  response;  the  filtered  signal  transform  is 
divided  by  the  instrument  transfer  function.  The  resulting  complex  spectrum 
is  then  Inverse  Fourier  transformed  into  the  time  domain,  to  produce  what 
will  hereafter  be  referred  to  as  the  filtered  signal. 

The  narrow  band  filtered  signal  will  appear  as  a quasi-sinusoidal 
carrier  wave  contained  within  a smooth  envelope.  The  next  step  in  the  program 
(Figures  14  and  15  ) is  to  construct  the  envelope  function  by  means  of  the 
Hilbert  transform.  In  particular,  a quadrature  signal  is  formed  by  multiplying 
the  transform  of  the  filtered  signal  by  -1  sgn(a))  and  then  Inverting  this  to 
the  time  domain  by  an  inverse  transformation.  The  envelope  function  is  then 
constructed  by  taking  the  square  root  of  the  sum  of  the  squares  of  the  filtered 
signal  and  its  quadrature.  The  maximum  of  the  envelope  function  occurs  at 
the  time  of  arrival  of  energy  at  the  center  frequency  of  the  filter  and  the 
amplitude  of  the  maximum  is  proportional  to  the  spectral  amplitude  of  the 
filtered  signal  at  the  center  frequency  of  the  filter.  The  instantaneous 
frequency  and  phase  are  computed  from  the  filtered  signal  and  its  quadrature 
signal  as  well  and  are  stored  for  subsequent  use  in  polarization  filtering 
with  additional  components  of  ground  motion. 

The  narrow-band  filtering  procedure  can  be  performed  on  a particular 
component  seismogram  at  a number  of  different  frequencies  within  some  band  of 
interest.  Correlation  of  the  resulting  envelope  functions  indicate.*!  the 
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arrival  times  of  the  various  frequency  components. 

The  results  of  a close-in  multiple  explosion  experiment  can  be 
used  to  desmontrate  the  ability  of  simple  narrow  band  filtering  to  achieve 
a separation  of  distinct  signal  arrivals  that  overlap  In  time. 

Specifically  we  use  a vertical  component  recording  from  a surface 
explosion  (Figure  16a)  to  simulate  a multiple  explosion  scenario.  This 
scenario  is  constructed  to  consist  of  a linear  array  of  three  equivalent 
yield  events,  equally  spaced  and  detonated  simultaneously.  To  simulate 
this  explosion  array  we  simply  delay  and  sum  the  recording  In  Figure  16a. 

The  delays  are  based  on  the  spacings  between  shots,  the  propagation  velocity 
assumed  over  this  spacing  and  the  azimuth  to  the  recording  station  (here 
taken  to  be  in  line  with  the  explosion  array).  The  resulting  composite 
seismogram  is  shown  in  Figure  16b  and  consists  of  three  separate  overlapping 
events,  each  consisting  of  a complex  wave  train. 

Both  the  original  and  composite  time  series  shown  in  Figures  16a  and  b 
were  narrow  band  filtered  with  sixteen  different  filters  with  center 
frequencies  varying  from  25  cps  up  to  100  cps,  since  the  basic  data  were 
digitally  recorded  at  a rate  of  500  saraples/second. 

Using  the  Hilbert  transform  procedure  mentioned  previously,  envelope 
functions  were  constructed  for  each  of  the  sixteen  filtered  outputs.  The 
envelope  functions  at  several  frequencies  are  shown  in  Figure  16c  and  d for 
the  original  and  composite  signals,  respectively.  First  we  observe  from 
this  array  of  filter  envelopes  in  Figure  16c  that  the  original  explosion  has 
three  separate  arrivals,  or  phases,  making  up  the  signal  wave  train. 

Secondly,  from  Figure  16d  we  can  see  the  triplication  of  each  of  these 
arrivals  corresponding  to  the  three  o:s.])]osiond  in  the  multiple  event  scenario. 
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The  clear  separation  of  the  signals  Is  even  more  dramatic  In  Figure  16f,  where 
the  sum  of  the  envelopes  at  the  different  frequencies  Is  plotted  fot  the 
composite  signal,  A similar  sum  Is  plotted  In  Figure  16a  for  the  original 
signal  which  again  shows  the  original  three  phases  of  the  single  event.  By 
contrast,  referring  to  Figures  16a  and  16b  , It  Is  not  at  all  obvious  from 
the  original  time  series  where  the  various  arrivals  occur  or  that  16b  is  a 
multiple  event. 

The  time  separatil  ons  of  the  maximum  power  arrivals  between  1.0  and  1.4 
seconds  in  Figure  16f  correspond  as  well  as  can  be  determined  to  the  time 
delays  associated  with  the  equal  spacings  of  the  three  simulated  events.  In 
addition,  note  that  the  amplitude  of  the  composite  signal  in  Figure  16b  is 
more  than  twice  the  amplitude  of  the  original  signal  in  Figure  16a,  but  that 
each  peak  in  the  composite  envelope  sura  is  nearly  equa]  to  the  peak  for  the 
original  signal  (Figure  16e).  Hence  the  power  in  each  of  the  three  separate 
signals  can  be  resolved. 

The  power  of  narrow  band  filtering  procedures  for  seismic  applications 

is  also  apparent  in  problems  of  separating  seismic  signals  from  noise. 

Figure  17  shows  an  actual  recording  of  a small  explosion  signal  emersed  in 

noise,  shovm  in  the  first  figure  at  the  upper  left.  The  signal  to  noise 

ratio  is  very  low  and  the  presence  of  a signal  is  not  at  all  obvious.  The 

successive  sequence  of  time  series  shown  are  narrow  band  filter  outputs  for 

filters  at  the  various  frequencies  indicated  by  the  value  of  f^  shown  in 

the  insets.  Clearly  the  noise  at  f^  = .3  Hz  and  1.0  Hz  obscures  the  signal, 

while  at  f = 3.0  H , 4.5  H and  6.0  H the  signal  to  noise  ratio  is  near  10. 
c z z z 

This  shows  that  for  seismic  signals  and  noise,  narrow  band  filtering  at 
high  frequencies  can  be  very  effective  for  signal  detection,  isolation  and 
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Figure  17:  Examples  of  increasing  signal-to-noise  ratio  for  a presumed  explosion 
(top  left  frame)  recorded  at  Norway  with  successive  application  of  high  frequency 
narrow-band  filters.  Arrow  at  top  denotes  approxinuite  arrival  times  of  the  P-wave 
on  the  unfiltered  trace. 
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for  spectral  estimates  of  the  signal. 

In  the  event  that  multicomponent  ground  motion  data  are  available, 
polarization  filters,  designed  to  pass  ground  motion  having  the  appropriate 
particle  motion,  can  also  be  applied  to  the  filtered  data  for  the  purpose 
of  identifying  and  separating  waves  of  differing  polarization.  The  sensing 
of  the  polarization  of  the  wave  field,  in  combination  with  narrow  band 
filtering  is  an  appropriate  means  of  automatically  identifying  seismic 
phases  or  wave  types.  In  many  cases  when  the  signal  wave  train  is  a complex 
superposition  of  phases  emersed  in  seismic  background  noise,  this  approach 
is  practically  the  only  systematic  means  of  identifying,  isolating  and  timing 
seismic  phases. 


Figure  18  is  an  example  of  the  application  of  combined  narrow  band 
filtering  and  polarization  filtering  using  a three  component  recording  of 
an  underground  nuclear  explosion  (code  name  Chartreuse).  The  event  was 


recorded  at  Kanab,  Utah  (KNUT)  at  a distance  of  312  km.  Because  of  the  short 
distance  between  the  source  and  receiver  the  record  is  parti  .ularly  compacted 
and  complicated,  with  many  arrivals  compressed  into  a fairly  brief  time  interval. 
The  traces  (B) , (C)  and  (D)  are  examples  of  the  combined  narrow  band  filtering 
and  polarization  filtering  at  two  frequencies.  Compressional  (P)  and  vertically 
polarized  shear  (SV)  waves  are  shown,  as  well  as  Rayleigh  type  surface  waves 
(elllptically  polarized).  Horizontally  polarized  body  shear  waves  and  Love 
waves  can  also  be  extracted  using  both  horizontal  components  of  motion,  but 
are  not  shown  here.  The  arrival  times  predicted  on  the  basis  of  tabulations 
appropriate  to  the  path  (Gutenberg,  1944;  Evernden,  19f'7;  Herrin,  19f.8)  are 


indicated  on  the  filtered  traces.  There  it  should  be  noted  that  the  filtering 
corresponds  to  "zero  phase  shift"  fllterfng  ,,„d  consequently  the  body  phase 
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arrival  times  ("group  times")  are  at  the  maxima  of  the  envelopes  of  the 
filtered  traces,  not  at  the  onset  of  the  motion  of  each  Individual  pulse  on 
the  filtered  trace.  Comparison  of  the  filtered  trace  envelope  maxima  with 
the  predicted  body  wave  arrival  times  Indicates  very  close  agreement  with 
differences  being  of  the  order  of  .1  to  .7  seconds.  One  would  no*"  expect 
any  better  agreement  than  this  since  the  predicted  travel  times  are  averages 
rather  than  simple  path  predictions.  Some  of  the  other  arrivals  isolated  by 
the  filtering  are  multiply  reflected  phases,  such  as  PP,  PS  and  so  on,  while 
the  Rayleigh  waves  in  (D)  are  either  signal  generated  Rayleigh  wave  "noise" 
or  higher  mode  Rayleigh  waves  from  the  source. 

A principal  product  of  the  filtering  operations  embodied  in  this  signal 
analysis  program.  In  addition  to  the  signal  isolation-phase  detection  and 
identification  output,  are  the  estimates  of  event  magnitude.  These  routinely 
consist  of  the  ordinary  body  wave  magnitude,  m^,  and  surface  wave  magnitude 
Mg  measured  at  1 Hz  and  .05  Hz  respectively  and  the  variable  frequency 
magnitudes,  mj^(f),  Mg(f)  (also  denoted  m^  and  M^)  corresponding  to  body 
and  surface  wave  magnitudes  measured  at  a variety  of  frequencies.  These 
latter  are  essentially  amplitude  spectral  estimates,  at  a set  of  discreet 
frequencies,  of  particular  wave  types  where  the  amplitudes  have  been 
converted  to  a magnitude  index  by  simply  dividing  by  the  appropriate  period 
T (i.e.,  the  center  period  of  the  narrow  band  filter)  and  taking  the  base 
ten  logarithm  of  the  result.  In  addition,  a constant  distance  correction 
factor  is  often  also  applied  to  give  magnitude  indices  for  the  event  that  are 
comparable  with  the  usual  m^^  and  magnitudes. 

Thus  magnitudes  may  be  computed  at  a variety  of  frequencies  from  the 
narro-.-  hand  filti.i  output  and  r.Jnco  ve  r;ay  also  associate  r.  group  time  or 
energy  arrival  time  with  each  a::’plitudc  measurement,  we  can  identify  the  wave 
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type  with  which  the  magnitude  measurement  is  associated. 

Figure  19  is  an  example  of  the  application  of  the  nijj(f)  data,  generated 
from  a large  number  of  events,  to  the  discrimination  of  earthquakes  and 
explosions.  The  event  data  shown  are  only  two  of  the  frequency  dependent 
magnitudes  that  characterize  the  events,  but  for  the  teleseismlc  distances 
(A  = 60-70°)  and  the  type  of  receiver  (LASA  seismic  array)  the  two  frequency 
dependent  magnitudes  at  .45  and  2.25  Hz  provide  the  best  separation  of  the 
explosion-earthquake  populations  in  the  ni|j(f)  parameter  space. 

The  large  scatter  of  the  earthquake  population  in  the  ni^(f)  space  is 
due  to  variability  of  the  stress  drop,  rupture  dimensions  and  rupture  rate 
among  the  events.  The  earthquakes  nearest  the  explosion  population  have  the 
highest  stress  drop  and  rupture  rate.  The  events  with  the  largest  values 
have  the  largest  failure  zone  dimensions. 

Some  of  the  scatter  in  the  event  populations  for  both  the  earthquakes 
and  explosions  is  due  to  background  seismic  noise.  This,  of  course,  is 
relatively  more  important  for  events  of  small  magnitude,  since  the  signal 
power  will  be  low  relative  to  the  noise  povjer,  and  causes  some  mixing  of  the 
populations  at  the  very  low  magnitudes.  (The  smallest  explosions  shown  have 
conventional  magnitudes,  m^^,  of  about  4). 

In  view  of  the  desirability  of  estimating  noise  properties  and  then 
correcting  for  noise  contamination,  the  signal  analysis  program  normally 
obtains  noise  spectral  estimates  and  uses  these  to  obtain  unbiased  estimates 
of  the  signal  spectral  magnitudes.  (In  this  case,  the  average  noise  spectral 
amplitude  at  the  appropriate  frequency  is  subtracted  from  the  signal  spectral 
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amplitude.)  The  nature  of  the  noise  has  been  investigated  for  each  of  the 
events  making  up  the  event  set  shown  in  Figure  19  using  the  signal  analysis 
program  and  it  was  found  that  noise  time  series,  when  viewed  in  an 
plane  such  as  that  in  Figure  19,  occupies  the  same  region  as  do  the  small 
magnitude  earthquakes.  This  means  that  when  the  noise  contaminates  a small 
explosion  signal,  the  result  will  be  such  as  to  make  the  explosion  appear  to  be 
more  ear thquake- like  and  for  low  signal  to  noise  ratios  to  actually  move  the 
explosion  point  into, or  very  near, the  earthquake  population.  Further  the 
noise  contamination  will  move  the  earthquake  points,  but  always  within  the 
earthquake  population  and  not  outside  of  it. 

This  is  illustrated  in  Figure  20  where  a number  of  the  small  events 
observed  at  LASA  and  shown  in  Figure  19  have  been  detected  at  a Norway 
array  and  processed  and  plotted  in  the  plane  for  f = . 6 Hz  and  f = 5.  Hz. 

Only  the  smaller  magnitude  events  are  shown  in  Figure  20.  Again,  in  this  plane, 
we  observe  convergence  of  the  populations.  The  arrows  on  the  smallest  explosion 
data  points  indicate  the  direction  in  which  a noise  correction  will  move  these 
points  in  this  magnitude  plane.  Figure  21  shows  the  effect  of  the  noise 
correction  applied  to  all  the  events,  including  earthquakes  as  well  as 
explosions.  (In  this  figure  only  the  log  of  the  amplitude  is  plotted,  the 
m^  plot  would  be  identical  except  for  a scale  change.)  It  is  evident  that 
the  explosion  population  is  now  well  separated  from  the  earthquake  population. 

The  arrows  on  the  two  small  explosions,  both  with  regular  m^  magnitudes 
somewhat  less  than  A.O,  indicate  that  it  is  likely  that  the  unbiased  estimate 
of  the  noise  contamination  is  lower  than  that  actually  present  and  that  the 
true  event  points  are  still  lower  in  the  direction  indicated. 
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Figure  20:  Spectral  magnitude  estimates  at  = 0.6  Hz  and  = 5.0  Hz  for  an  event 
population  recorded  at  the  Oyer  array  in  Norway. 
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In  any  case,  it  is  evident  that  the  combined  signal  and  noise  analysis 
capability  is  important  for  the  correction  of  signal  data  in  order  to  obtain 
better  estimates  of  source  characteristics,  as  expressed  by  the  observed 
signal  information. 

In  total,  the  current  signal  analysis  program  constitutes  an  extremely 
useful  and  quite  powerful  tool  to  be  used  for  the  detection,  identification 
and  analysis  of  seismic  signals.  The  output  of  this  processor  provides  not 
only  the  required  data  for  explosion-earthquake  discrimination  but  also 
quantitative  noise  corrected  data  for  analytically  based  determinations 
of  source  and  medium  properties,  such  as  stress  drop  and  rupture  dimensions 
of  earthquakes,  energy  yields  for  explosions  and  travel  times  and  spectral 
amplitudes  of  body  and  surface  waves  for  medium  studies.  Thus  it  constitutes 
a low  cost,  rapid,  automated  procedure  for  the  detailed  analysis  of  very 
large  event  data  sets. 


V.  Stress  Estimatas  Along  the  North  Pacific  Circumference. 


Observational  data  has  been  used  to  infer  stre:-.  • levels 

in  the  Alaskan-Aleutlan  and  Japan-Kuri  les-Kamcliatka  rep.ions  u'iint'  • 
approach  discussed  in  the  previous  sections.  Figure  22 
shows  all  the  data  available  for  the  Alaskan  region,  plotted  against 
theoretical  curves  for  dip-slip  events  in  the  0 - 20  kni  depth  range. 

Only  some  of  this  shallow  event  data  can  be  interpreted  in  terns  of  the 
theoretical  curves  shown  of  course,  but  the  curves  do  indicate  the 
stress  range  for  the  variety  of  events  in  the  depth  range  from  0 to  50  km. 

Similar  data  for  the  Calif ornia-Nevada  region  is  sho^vm  in  Figure 
(23).  Clearly,  the  amount  of  available  data  for  this  region  is  less 
than  that  for  the  very  active  Alaska-Alcutian  region.  Consequently, 
for  stress  estimates,  the  present  data  that  v.’e  have  must  be  supple- 
mented by  additional  observations. 

Figure  (24)  shows  the  currently  available  data  for  the  Japan  region, 

plotted  in  the  same  manner  as  the  data  for  the  other  regions.  This 

region  is  the  most  active  in  the  world  and  provides  an  indication  of 

the  maximum  amount  of  m,  - M data  that  can  be  obtained  over  a five 

b s 

to  six  year  period  from  the  U.S.G.S.  compilations.  The  data  for  the 
Japan  and  Alaskan  regions  are  sufficient  for  an  application  of  the 
approach  outlined  in  our  earlier  discussions. 

At  this  point  it  is  evident  that  we  can  obtain  stress  drops  and 
fault  dimensions  by  direct  comparison  of  theoretical  predictions  with 
the  data,  as  indicated  in  Figures  (22)  to  (24).  The  data  must  however 
be  separated  into  groups,  corresponding  to  activity  in  particular  depth 
zones  for  earthquakes  of  a particular  type,  with  these  event  subsets  inter- 
preted using  comparable  theoretical  results. 


SURFACE  WAVE  MAGNITUDE 
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Figure  22:  Observed  m,-M,  data  for  the  A1  ska-Alcutian  Arc  Seismic  Region  shc-m  with 
theoretical  m, -M , curves  with  prestress  as  a parameter  for  the  curves  and  fault  length 
variable  along  each  curve.  Theoretical  curves  and  data  are  for  dip  slip  and  thrust 
events . 
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Figure  23:  Observed  . data  for  the  California  Nevada  seismic  region  shown  with 

theoretical  m,-M  curves  ?or  thrust  and  dip  slip  type  earthquakes. 
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Figure  24:  Observed  data  for  the  Japan-Kur  i 1 es-Kamelia  t ka  seismie  region  sliown 

with  theoretical  m,-M^  curves  for  thrust  or  dip  slip  earthquakes.  The  events  observed 
are  essentially  all  thrust  or  dip  slip  type.  The  mean  stress  drop  fi^r  tlie  population  is 
somewhat  larger  than  100  bars  within  any  events  having  stress  drops  near  1 kilohar. 
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if  we  assume  that  melting  occurs  during  or  as  a consequence  of 
failure  within  the  rupture  zone,  then  the  stress  drop  is  in  fact  equivalent 
to  the  recoverable  ambient  stress  level.  We  need  not  make  this  assumption 
however  in  that  we  can  use  the  largest  stress  drop  events  to  provide 
the  estimate  of  the  ambient  stress,  under  the  weaker  assumption  that 
at  least  some  failure,  processes  lead  to  partial  melting  (or  similar  process 
leading  to  a total  lack  of  shear  strength  ) and  hence  a "total  stress 
drop"  on  the  failure  boundary.  We  can  then  use  the  event  stress  drop 
estimates  from  a large  population  of  events  to  infer  the  ambient 

stress  by  plotting  all  the  events  in  a restricted  depth  range  on  a map; 
and  to  then  contour  the  stress  sampling  using  the  largest  stress 

drop  events  to  control  the  contours.  Defining  the  ambient  stress  estimate 
in  this  way  then  allows  us  to  determine  whether  the  estimated  stress 
from  all  the  events  conforms  to  the  stress  level  pattern  defined  by 
the  highest  stress  drop  events.  If  all  the  stress  drop  estimates  do  in 
fact  consistently  conform  in  value  with  the  contour  levels,  then  we  could  reason- 
ably infer  that  essentially  all  the  events  had  stress  drops  from  the. 
ambient  level  to  the  same  fixed  level  and  that  most  likely  this 
stress  drop  was  total,  with  partial  melting,  or  a similar  weakening 
process,  involved.  This  conclusion  is  of  course  predicated  on  there 
being  a large  data  sample — a small  sample  would  obviously  disallow  such 
a conclusion.  In  addition,  it  is  required  that  no  very  large  event  take 
place  in  the  time  window  of  the  data  sample,  since  it  would  so  strongly 
perturb  the  stress  as  to  cause  subsequent  events  in  the  vicinity  to  have 
much  different  stress  drops. 
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Results  from  the  A 1 askauj^ leutiau  Rep.lon 

Figures  (Z5  ) through  (28)  show  the  results  ol  applying  tliis  procedure 
to  a fairly  large  total  population  of  events  from  the  Alaskan-Aleutian 
region.  Hero  the  events  were  separated  into  groups  based  on  depth  and 
type  so  that  the  stress  estimates  are  plotted  in  particular  depth  ranges. 

These  figures,  taken  together,  present  a picture  of  the  stress  level  as 
a function  of  depth  for  the  epicentral  region  covered  by  the  zone  outlined 
as  "Seismic  Region  1".  The  contouring  was  made  using  the  highest  stress 
drop  events,  as  previously  described.  A clo.se  look  at  the  individual 
stress  drops  for  the  whole  population,  (see  the  figure  legend  for  an  explana- 
tion of  the  plotted  numbers  associated  with  each  event)  relative  to  the 
contours  shows  that  with  very  few  exceptions  (i.e. . two  or  three  events) 
all  the  events  are  consistent  with  the  contoured  stress  levels.  This 
suggests  (but  does  not  necessarily  require)  that  all  the  stress  drops  reduce 
the  stress  to  a fixed  level  at  the  failure  boundary  or  to  a fixed  fraction 
of  the  ambient  level.  In  view  of  the  wide  variations  in  material  properties 
to  be  expected  both  in  depth  and  laterally  along  the  plate  boundary 
however,  such  consistency  would  strongly  imply  a total  stress  drop  for  the  events. 

Certainly  one  of  the  striking  features  of  these  results  is  the 
consistent  variation  of  the  stress  with  depth  and  along  the  plate  boundary. 

In  addition,  the  data  give  very  high  stresses  in  the  region  from  Kodiak 
Island  to  the  end  of  the  Alaskan  peninsula  and  this  picture  is  consistent  with 
the  presence  of  a seismic  gap  in  this  region  (eg^.  Sykes,  1971)  and  the  fact 
that  the  recurrence  data  for  large  earthquakes  in  this  region  implies  another 
large  earthquake  zone  in  the  near  future  (ej^._  Sh.alcal  and  Willis,  1972). 

Further  we  find  that  this  stress  estimation  methods  gives  rather  low  stress 
levels 
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for  the  zone  of  the  1964  Alaskan  event,  wliich  would  be  expected. 
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Figure  25:  Epicenter  stress  map  from  events  in  the  Alaskan-Al eu t ian  Arc  region  having 
focal  depths  in  the  range  from  0 to  20  km.  Each  event  is  denoted  by  a time  sequenced 
event  number,  the  stress  drop  and  its  m^^  and  values  as  shown  in  the  legend.  The  stress 
drop  values  are  contoured  to  show  the  spatial  stress  pattern  and  provide  a useful  indica- 
tion of  the  consistency  of  the  results.  The  event  numbers  provide  an  indication  of  the 
time  sequence  of  the  sampling.  If  shear  molting  occurs  in  the  failure  zone  (on  the 
failure  "plane"),  then  the  stress  levels  shown  are  equal  to  the  ambient  prestress  levels 
in  the  material. 
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Figure  2h:  Epirt-ncral  stress  m;ip  from  events  in  the  A1  nsknn-A] eu t ian  Arc  region  having 
focal  depths  in  the  20  to  33  km  range.  Events  listed  as  being  at  a depth  of  33  km  are 
omitted  from  the  data  set  shown  and  are  treated  separately  since  the  assignment  of  a 33 
km  depth  is  the  indication  of  unknown  depth  for  an  event.  Note  the  high  stress  zones 
along  the  Alaskan  peninsula  in  this  depth  range.  This  high  stress  zone  appears  to  be 
correlated  with  the  high  stress  region  plunging  downward  and  eastward  along  the  Pacific 
plate  houndary. 
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Figure  27:  Eplcentral  stress  map  from  events  in  the  Alaskan-Aleutian  Arc  region  having 
focal  depths  in  the  33  to  40  km  range.  Events  listed  with  hypocenters  at  33  km  are 
omitted.  The  high  stress  zone  along  the  Alaska  Peninsula  is  seen  to  persist  through 
this  depth  interval,  with  the  stress  levels  being,  again,  near  the  1 kilobar  level. 
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Figure  28:  Epicentral  stress  map  from  events  in  the  Alaskan-Aleutian  Arc  region  having 
focal  depths  in  the  40  to  50  km  range.  The  high  stress  zone  along  the  peninsula  is 
larger  and  more  extensive  in  this  depth  range  than  at  the  shallower  depths.  It  is  likely 
to  be  the  depth  rang"  for  nucleation  of  a very  large  earthquake  as  a result  of  failure 
along  the  plate  boundary  from  just  west  of  Kokiak  Island  to  the  end  of  the  Alaskan 
Peninsula.  For  stress  levels  near  1 kbar,  as  shown,  this  event  could  occur  at  any  time. 
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Results  from  the  Japan-Kur LlGS-Kainchatka  ReRion 

Figures  (29)  through  (3-2)  show  the  inferred  stress  levels  for  the 
Japan  region  in  the  depth  range  from  0 to  50  km.  Somewhat  more  data  is 
available  for  this  region  than  for  the  Alaskan  zone  and  the  stress  estimates 
are  correspondingly  somewhat  more  detailed — and  reliable.  An  outstanding 
feature  of  these  results  are  delineated  zones  of  very  high  stress,  in 
the  depth  range  20  - 50  km, of f the  coasts  of  Kamchatka,  Hokkaido  and  one 
off  Honshu,  near  Tokyo.  These  patterns  show  some  irregularity  with  depth 
but  in  general  the  high  stress  patterns  are  found  to  persist  over  a 
considerable  range  in  depth.  As  with  Alaska,  the  high  stress  zones 
are  generally  at  depths  greater  than  20  km,  with  the  stress  levels  near 
the  surface  more  moderate.  The  one  exception  appears  to  be  the  very 
high  stresses  indicated  in  the  Kuriles  region  in  Figure  (29).  This 
stress  zone  appears  to  be  located  in  the  North  American  plate,  north- 
west of  the  trench  marking  the  plate  boundary,  rather  than  entirely  in 
the  narrow  zone  along  the  boundary  itself. 

A particularly  critical  test  of  the  entire  approach  is  available 
using  the  data  set  for  this  region.  Specifically,  we  should  obtain  a 
high  stress  level  in  a concentrated  zone  before  a large  earthquake,  while 
after  the  event  we  should  obtain  a stress  pattern  which  is  considerably 
changed.  One  would  expect  a general  reduction  in  the  stress  with  perhaps 
some  high  stress  concentration  levels  at  the  ends  of  the  failure  zone. 

The  Japan  region  data  contained  two  large  events  with  reasonably  large 
numbers  of  small  events  occurring  both  before  and  after  them,  and  conse- 
quently we  can  test  the  seismic  stress  estimates  obtained  against  the 


criteria  that  the  local  stress  level  be  generally  lower  after  such  a 


large  event. 
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Figure  29:  Epicentral  stress  map  from  events  in  the  Japan-Kuriles-Kamchatka  region  for 
the  depth  range  from  0 to  20  km.  The  high  stress  zone  in  the  Kuriles  region  corresponds 
to  a shallow  stress  concentration  of  apparently  large  dimension  within  the  Asian  plate, 
which  is  being  underthrust  by  the  Pacific  plate. 


Figure  30:  Epicentral  stress  map  from  events  in  the  Japan-Kuriles-Kamchatka  region  for 
the  depth  range  from  20  to  33  km.  High  stress  zones  occur  along  the  Kamchatka  and 
Hokkaido- Honshu  coasts,  with  the  latter  zone  showing  a complex  pattern  of  high  and  low 
stresses  characteristic  of  zones  of  recent  large  earthquakes  such  as  event  no.  61. 
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Figure  31:  Epicentral  stress  map  from  events  In  the  Japan-Kuriles-Kamchatka  region  for 
the  depth  range  from  33  to  40  km.  The  high  stress  zone  off  the  Honshu  coast  appears  more 
regular  and  extensive  In  this  depth  range  and  constitutes  a high  risk  zone-especlally 
In  view  of  Its  proximity  to  Tokyo. 
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Figure  32:  Epicentral  stress  map  for  the  Japan-Kuriles-Kamchatka  seismic  region  from 
events  having  focal  depths  in  the  depth  range  40  to  50  km,  with  events  occurring 
before  the  large  event,  No.  193,  represented.  This  event  and  its  aftershocks  occurred 
within  the  moderately  high  stress  zone  off  the  northeast  point  of  Hokkaido  and  results 
in  a considerable  change  in  the  stress  level  and  spatial  pattern. 


72 


Figure  (33)  shows  the  stress  magnitude  distribution  for  the  40  to 
50  km  depth  range  after  one  of  the  large  events  (No.  193)  that  would 
be  expected  to  strongly  perturb  the  stress  field  in  a large  zone.  We 
would  expect  to  observe  a considerable  change  in  the  stress  field 
as  obtained  from  the  sampling  provided  by  the  subsequent  small  events. 
Comparison  of  Figures  (32)  and  (33)  shows  a major  shift  in  the  stress 
pattern  due  to  the  large  event,  with  the  high  stress  zone  shifting 
landward  after  the  event  and  a considerable  reduction  in  stress  in  the 
immediate  vicinity  of  the  large  earthquake.  Hence,  our  spatial  resolution 
of  the  stress  appears  just  sufficient  in  this  case  to  verify  that  the 
stress  change  is  as  would  be  expected. 

Figures  (34)  and  (35)  show  the  stress  field  contours  and  the  event 
data  before  and  after  another  large  event  (Event  No.  144).  This  event 
occurred  just  at  the  point  where  the  Aleutian  arc  intersects  the  Kamchatka 
Peninsula.  Figure  (34)  shows  that  prior  to  the  event  the  zone  was  rather 
highly  stressed,  with  the  stress  level  reaching  at  least  500  bars. 

After  the  event.  Figure  (33)  shows  that  the  later  small  events  provide 
a sampling  of  a stress  field  that  is  generally  very  much  lower  in  this 
region,  with  however,  some  very  small  zones — presumably  near  the  ends 
of  the  newly  failed  zone  associated  with  the  large  earthquake — where 
the  stress  is  quite  high.  Thus  the  time-space  picture  presented  by  this 
event  sequence  again  entirely  conforms  to  what  we  know  should  occur. 

Finally  we  also  observe,  as  for  the  Alaska  region,  that  the  stress 
levels  obtained  by  contouring  on  the  basis  of  the  largest  stress  drop 
values  obtained  at  each  point  are  consistent  with  (very  nearly)  all  the 

j, 

f observations.  That  is,  when  a high  stress  zone  is  delineated  by  a 
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Figure  34:  Epicentral  stress  map  for  the  Japan-Kur iles-Kamchatka  seismic  region  from 
events  having  focal  depths  given  as  33  km,  with  only  the  events  before  the  large  event 
No.  144  represented.  Event  No.  144,  which  occurred  off  the  coast  of  Kamchatka,  and  its 
aftershocks  are  not  included  so  that  the  stress  field  prior  to  this  event  is  shown. 

It  corresponds  to  the  northern-most  high  stress  zone  along  the  Kamchatka  coast. 
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Sampling  of  events  for  a region,  then  we  find  that  subsequent  (small) 
events  in  the  zone  will  all  be  high  stress  drop  events  until  a large 
event  occurs,  at  which  time  the  stress  field  is  greatly  changed,  usually 
over  the  entire  zone. 

We  conclude  that  the  general  space-time  picture  obtained  from  this 
stress  estimation  procedure  is  quite  consistent  with  the  expected  stress 
behavior  within  the  active  seismic  zones  studied.  Further,  the  results 
strongly  suggest  that  most  events  involve  total  (shear)  stress  drops 
on  the  failure  boundary  so  that  the  stress  drop  estimates  are  a direct 
measure  of  the  recoverable  nonhydrostatic  stress  field  in  the  medium. 
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VI.  Summary  and  Conclusions. 

Some  conclusions  that  have  been  drawn  that  are  oE  importance  for  the 
earthquake-explosion  discrimination  problem  and  for  a p.cneral  understanding 
of  earthquake  elastic  v;ave  radiation  fields  are:  ; 

(i)  A relaxation  source  theory  appears  to  provide  predictions  of  the  i 

! 

radiation  field  from  earthquakes  that  are  in  good  agreement  v/ith  obser- 
vations,  specifically  and  other  magnitude  data  and  spectra  from 

near  field  observations  of  earthquake  radiation. 

(ii)  The  earthquake  theory,  combined  with  explosion  source  model  theory, 

provides  a coherent  explanation  of  the  discrirainent.  It  also 

provides  a predictive  basis  for  ascertaining  v/hen  and  under  what  cir- 
cumstances this  discriminent  can  fail.  That  is,  it  provides  a basis  of 
understanding  of  anomalous  events. 

(iii)  Theoretical  predictions  of  spectral  magnitude  discriminents 
indicate  that  a short  period  body  wave  magnitude  discrimination  can  be 
at  least  as  effective  as  discrimination  and  in  general  provides  a 

wider  separation  of  explosion-earthquake  populations  to  lower  magnitudes 
(m^  V 4.0)  vzith  few  if  any  anomalous  shallow  earthquakes  appearing  in 
the  explosion  population.  Observational  tests  of  the  theoretical  pre- 
dictions using  roughly  200  Eurasian  events  have  confirmed  these  predictions. 

The  only  earthquakes  within  the  explosion  population  were  some  very  deep 
earthquakes.  In  addition,  both  single  and  multiple  explosions  (the  latter 
designed  t appear  as  earthquakes  to  all  other  discriminents,  including 
vjere  identified  by  this  discriminent  as  explosions.  This  dis- 
criminate is  easily  Implemented  in  the  field  and  as  an  on-line  processor 
to  detect  and  identify  event  signals. 
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Other  spectral  magnitude  measurements  employing  surface  waves  were 
also  shown  to  be  discriminatory,  but  the  theoretical  predictions  have  not 
been  observationally  verified  in  detail. 

(iv)  Analysis  of  several  hundred  events  in  the  North  Pacific  region  shows 
that  the  stress  environment  surrounding  a failure  zone  and  leading  to 
failure  can  in  many  cases  be  highly  inhomogeneous  with  local  nonhydro- 
static stress  concentrations,  of  relatively  small  spatial  extent,  of  the 
order  of  1.0  kb.  Most  moderate  sized  earthquakes,  however,  appear  to 


occur  in  zones  with  average  prestress  at  the  100  bar  level,  with  the  stress 
quite  uniform  over  regions  of  the  order  of  100  km  in  radius. 


(v)  Near  field  spectra  observations  and  data  for  small  earthquakes 


(nij^  ^ 5),  suggest  strongly  peaked  far  field  spectra  for  at  least  some 
earthquakes.  A conclusion  to  this  effect  is  compatible  with  the  deduced 
strong  spatial  variations  in  stress  (indeed  it  would  logically  follow) 
provided  that  the  Rg  factor  approximation  used  in  the  relaxation  source 
theory  is  a reasonably  accurate  meai.s  of  representing  the  radiation  from 
relaxation  of  spatially  variable  prestress.  A consequence  of  strongly 
peaked  far  field  radiation  is  a reduction  in  the  observed  value,  with 
the  event  appearing  explosion-like  in  the  plane. 

(vi)  Comparisons  of  theoretical  predictions  from  the  analytical 
relaxation  source  theory  models  and  complex  purely  numerical,  two  and 
three  dimensional  earthquake  models  employing  failure  criteria,  including 
plasticity,  are  in  substantial  agreement.  This  agreement  lends  strong 
support  to  the  reliability  of  the  "first  order"  relaxation  theory 
models  of  earthquakes. 

(vil)  A general  Green's  function  representation  of  the  radiation 
field  produced  by  rapid  material  failure  in  a generally  prestressed 
medium  has  been  obtained.  The  theory  includes  the  complex  dynamical 
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boundary  conditions  for  an  expanding  failure.  It  is  shown  that  the 
theory  is  amenable  to  approximations  allowing  reliable  and  quite 
c general  models  to  be  constructured  and  for  which  analytical  results 

can  be  obtained.  Development  of  these  higher  order  models  was 
initiated  during  the  period  of  this  study  with  some  preliminary  results 
obtained  by  the  end  of  the  contract  period. 
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Abs^tjict 

Spontaneous  failure  In  n solid  medium  is  described  as  a localized 
transition  of  the  material  from  one  physical  state  to  another, 
characterized  in  part  by  contrasting  rheological  properties  and 
density.  Such  a process  is  viewed  as  a local  disordering  of  the 
relatively  ordered  structure  of  the  solid  due  to  any  variety  of  causes, 
such  as  massive  microfracturing  or  shear  melting,  and  can  be  confined  to 
a very  thin  zone,  but  nevertheless  of  finite  volume  such  that  a volumetric 
transition  energy  can  be  defined.  This  leads  to  the  description  of 
failure  as  a generalized  phase  transition  in  a prestressed  continuum, 
with  instability  and  transition  zone  growth  being  driven  by  the  energy 
contributions  from  the  relaxation  of  stress  in  the  surrounding  medium. 

Direct  application  of  mass,  momentum  and  energy  conservation  to  such  a 
generalized  phase  transition  leads  to  "jump"  conditions  specified  on  the 
growing  boundary  surface  of  the  transition  zone,  that  relate  the  rupture 
growth  to  discontinuous  changes  in  the  dynamic  field  variables  across 
the  failure  zone  boundary.  These  field  discontinuities  are,  in  turn, 
related  to  the  localized  changes  in  physical  properties  induced  by 
failure.  Dynamical  conditions  for  rapid  spontaneous  failure  growth  in 
a stressed  medium  arc  investigated  in  some  detail,  and  we  find  that  the 
failure  boundary  growth  can  be  simply  expressed  in  terms  of  energy 
"failure  condition"  and  a dynamic  growth  condition  specifying  the 
rupture  velocity.  These  results  imply  that  the  internal  energy  change 
associated  with  earthquakes  is  in  the  range  10^-10^  ergs/gm.  Further 
the  failure  growth  rate  is  shown  to  be  expressible  in  terms  of  the 
rheological  properties  of  the  material  before  and  after  failure.  For 
shear  melting  resulting  in  a low  viscosity  fluid,  for  example,  the 
rupture  velocity  will  be  near  the  shear  velocity  of  the  original 
material.  A general  Green't  function  solution  for  the  radiation  due 
to  stress  relaxation  in  the  medium  surrounding  the  growing  failure 
zone  is  given  and  provides  the  basis  for  detailed  computations  of  the 
strain  or  displacement  field  changes  due  to  spontaneous  failure  j 

processes.  In  particular,  it  is  shown  that  the  jump  conditions  for  the 
growing  transition  zone  boundary  appear  naturally  as  surface  integral 
terms  over  the  boundary.  Since  these  boundary  conditions  contain  the 
failure  rate  explicitly,  then  these  terms  Include  effects  that  have 
not  been  represented  in  previous  integral  representations  of  the 
radiation  field  resulting  from  failure.  Further,  we  show  that  the 
formal  Green's  integral  representation  for  the  dynamical  wave  field  can 
be  used  with  known,  simple  Green's  functions  to  generate  approximate 
solutions  for  complex  failure  processes  occurring  in  media 
with  inhomogeneous  material  properties  and  prestress. 
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1.  Introduction 

A general  description  of  material  failure  under  stress  loading  is 
obtained  witliin  the  frame  of  continuum  mechanics  by  considering  the 
failure  process  to  be  a generalized  phase  transition  (Archambeau,  1969). 
In  this  description  the  failure  zone  occupies  a finite  volume,  the 
boundary  of  which  corresponds  to  a moving  phase  boundary  in  a 
prestressed  medium.  The  growth  of  the  failure  zone  is  then  controlled 
by  the  energy  balance  at  the  "phase  transition".  That  is,  we  can 


define,  as  a macroscopic  equivalent,  a "latent  heat"  of  transition, 
associated  with  the  energy  required  for  microfracture  and/or  partial 
melting  along  grain  boundaries  or  fracture  surfaces.  This  energy 
term  can  be  used  as  an  intrinsic  characteristic  of  the  material, 
and  describes  the  energetics  of  the  irreversible  processe;  of 
failure . 

From  this  point  of  view  it  is  not  necessary  to  specify  the 
microscopic  details  of  the  failure  process,  which  may  involve 
brittle  fracture,  plastic  flow,  melting,  etc.  ...,  but  only  the 
composite  irreversible  energy  required  for  these  processes,  as  a 
function  of  temperature,  pressure  and  material  type. 

The  dynamical  evolution  of  the  failure  zone  is  then  controlled 
by  the  laws  of  conservation  of  mass,  momentum,  and  energy,  in  a 
continuum  including  a moving  phase  boundary.  Because  of  the 
discontinuous  behavior  of  the  medium  at  the  failure  boundary,  the 
momentum  and  energy  equations  are  coupled  through  boundary 
conditions  at  the  failure  surface.  However,  one  can  still  construct 
a Green's  tensor  representation  theorem  for  the  radiated  field,  in 
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wliich  the  coupling  appears  as  a term  involving  rupture  velocity  in 
the  boundary  surface  integral  over  the  failure  surface.  Since  the 
radiation  field  arises  from  the  relaxation  of  stress  in  the 
surrounding  medium,  the  elastic  field  representation  corresponds 
to  that  for  relaxation  source  theory  (e.g.,  Archambeau,  1968,  1972; 

Minster,  1973). 

In  this  paper  we  shall  give  a concise  account  of  the  theoretical 
development  of  the  basic  conservation  conditions  that  must  hold  in  a 
continuum  with  localized  discontinuous  properties.  We  shall  then 
focus  on  the  case  where  the  dynamical  deformation  of  the  surrounding 
medium  induces  rapid  growth  of  a failure  zone,  since  it  corresponds 
to  observed  tectonic  failure  resulting  in  earthquakes.  (The  theory 
is,  of  course,  also  applicable  to  a variety  of  pliyslcal  processes, 
including  shock  waves  in  solids  and  ordinary  phase  changes.)  In 
this  regard  we  will  investigate  the  constraints  placed  on  the  rate 
of  failure  by  the  conservation  relations  ("jump  conditions")  that 

must  hold  on  the  rupture  boundary.  i 

Finally,  the  representation  of  the  radiation  field  which  arises  i 

from  stress  relaxation  around  the  failure  zone  is  expressed  in  terms 

1 

of  a general  Green's  function  solution.  In  order  to  obtain  such  a 
solution,  the  usual  linearizing  approximations  are  needed  (e.g., 
the  assumption  of  infinitesimal  strains  in  the  medium  exterior  to 
the  failure  zone),  as  well  as  some  approximations  that  are  unique 
to  this  problem.  We  will,  however,  develop  an  integral  representation 
of  the  radiated  field  with  the  minimum  number  of  assumptions  in 
order  to  obtain  general  results. 


I 
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2.  Transport  theorems  and  conservation  equationH 

In  a general  continuum  representation  of  material  flow  the 
conservation  of  a function  of  the  flow  F may  be  expressed  as 

^ — J F(x,t)  d X - J k(x,t)  d X , 


I'd) 


Uft) 


where  the  source  density  has  the  general  form 

k(x,t)  - K(x,t)  + V*  J(x,t) 

with  ^ the  flux  of  Interacting  fields  and  K corresponding  to 
the  Intrinsic  production  of  F In  the  volume  l'(t) 

Let  S(t)  be  the  boundary  of  the  material  volume  l/ft)  , 
let  material  velocity  field,  and  assume  that 

l/(t)  is  cut  by  a surface  of  discontinuity  I(t)  ; moving  with 
the  velocity  U(2i>t)  ^ as  Illustrated  In  Figure  1.  Let 
denote  the  jump  of  a quantity  G across  I . Then  the  Reynolds 
transport  theorem  takes  the  form 


^ f 

(/(t) 


t)  d^x 


I 

i'(t)0{i/(t)ni:(t)} 


+ d X 


4 


f Cf(  U-V  )•  nj 

i'(t)ns(t) 


da 


(1) 


(2) 


(3) 


Figure  1.  Case  of  a propagating  discontinuity.  r(t)  is  the 
surface  of  discontinuity  moving  with  velocity  U . 

S, (t)  + S,(t)  is  a material  surface,  moving  with  the  medium 
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and  Gauss's  theorem  takes  the  form 

^ V*^  d^x  + <!*>  " ^ 

W0{zfY/}  zny  sQisnz} 

where  © and  0 denote  the  set  theoretic  difference  and  Intersection 
respectively.  A simplified  proof  of  these  results  can  be  found,  for 
example,  in  Eringen  (1973).  Edelen  (1962)  states  a similar  result 
to  (A),  but  neglects  to  exclude  the  points  of  E from  the  volume 
Integral,  so  that  the  term  in  ^'n  is  effectively  included  twice 
in  his  result.  Note  that  (3)  and  (4)  hold  only  as  long  as  none  of 
the  functions  f » v , ^ has  a singularity  on  E stronger  than 
a mere  discontinuity  (Minster,  1973).  A detailed  proof  of  these 
results  can  be  obtained  using  the  theory  of  distributions. 

If  generalized  functions  are  used, 

it  is  easy  to  see  that  (4)  may,  in  fact,  be  written  in  the  usual 
form 

^ 7*J  d^x  = ^ ^'n  da  . (5) 

V S 

Using  (3),  (4)  in  (1),  (2),  and  noting  that  (1)  holds  for 
any  volume  l'(t)  , we  have  away  from  discontinuities 


+ V^- (FV-J)  = K , (6) 
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which  is  the  differential  con^i^■rv.■lt  Ion  equation.  Since  this  has  to 
hold  at  all  points  away  from  T , (1)  reduces  to 

J"  [EF(y-U) ’n])  da  - ([j'njda  (7) 

vu)ni(t)  V(i)  i>.(t) 

so  that,  since  l^(t)  is  arbitrary,  we  must  have  at  all  points 
on  Z 


(8) 


which  is  the  boundary  condition  to  be  satisfied  on  E for 
conservation  of  F 

A result  of  similar  form  has  also  been  given  by  Freund  (1970) 
in  a somewhat  different  context.  In  a study  by  Snoke  (1976), 
motivated  by  results  equivalent  to  (8)  given  by  Minster  (1973), 
an  attempt  was  made  to  obtain  general  ".lump  conditions"  such  as 
(8)  by  a method  of  transformation  of  the  equation  of  motion  to  a 
moving  coordinate  frame.  The  approach  does  not  contain  the 

I 

essence  of  the  "moving  boundary  problem",  however,  namely  the  j 

discontinuous  behavior  of  the  field  variables,  and  yields  results 
which  are  not  conservation  relations  on  the  moving  boundary 

surface.  Burrldge  (1976),  however,  employs  boundary  conditions  ' 

that  are  approximately  the  same  as  those  given  here  for  the  special 
problem  he  considers. 


A 


r 


If  wo  doflno  V^*  » y^-U  ns  tlio  velocity  of  the  matorini, 

relative  to  tlie  boundary  Z , and  Introduce  tiu*  current  density 

of  F per  unit  time  per  unit  area  (flux  vector):  H ~ :l  • 

* * 

and  the  relative  flux  vector  through  V ~ J ; then 

the  conservation  laws  have  the  form 


ap 

at 


K 


(9) 


ILf*  n^]]^  - 0 


(10) 


where  K is,  as  before,  the  rate  of  production  of  F per  unit 
time  per  unit  volume. 

We  observe  that  if  the  boundary  Z is  a material  surface, 
or  if  y^  ‘n  = 0 , then  (10)  becomes 

([j*n]]j-  = 0 (11) 


which  is  the  usual  boundary  condition  when  the  boundaries  move 
with  material  flow  or  deformation. 

Conservation  equations  apply  in  general  to  tensor  components, 
and  It  is  straightforward  to  show  that  the  general  forms  for  (9) 
and  (10)  are 


ar.  , 

9t 


9 

9x 


[F, 


,3m 


] = K 


i. . 


. jm 


= 0 


(12) 
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Table  1 furnishes  a sununary  of  the  most  common  conservation  equations. 
Additional  conservation  equations  are  also  possible,  sucli  as  those 
investigated  by  Fletcher  (1974)  for  hyperclastic  media  using  Noether's 
theorem.  The  additional  relations  involve  internal  angular  momentum 
for  polar  media,  and  since  we  are  only  concerned  with  nonpolar  media, 
conservation  of  angular  momentum  only  requires  that  the  stress  tensor 
T^j  be  symmetric.  Thus,  equations  (12),  together  with  Table  1, 
provide  the  complete  set  of  conservation  relations  for  the  media 
of  importance  in  this  study. 

An  alternative  expression  of  the  conservation  equations  in 
four  dimensional  form  is  also  possible.  We  find  such  a compact 
representation  to  be  useful  when  we  consider  integral  representations 
of  the  linearized  equations  of  motion  for  tlie  continuum  in  a 
following  section.  Here  we  simply  display  the  forms  for  the 
general  (nonlinear)  case. 

In  particular,  the  conservation  relations  can  be  put  in  a 
compact  explicit  form  that  expresses  conservation  of  mass  and 
momentum  by  introducing  fields  and  Kg  with  the  Greek 

indices  running  over  the  range  1 to  4.  Thus,  with  space-time 
coordinates  represented  by  x^  we  have  for  the  momentum  equation 


Here  we  retain  a non-relatlvistic  description,  and  merely  consider 
time  t as  a fourth  coordinate  x^  in  the  Newtonian  sense.  The 
summation  convention  over  all  repeated  indices  applies  to  both 
Greek  and  Roman  indices. 


1 


¥ 


loo 


3F 


a6 


9x 

a 


u.B  » 1,2, 3, 4 


(13a) 


*'^6  "a  ij  ■ » 


(13b) 


where  Kg  and  are  the  "space-llke"  variables  given  by 


Kg  - (pbj,  pb^,  pb^,  0)  = (pbj  , 0) 


(13c) 


"a  ' ^”l’  ”2’  "3’ 


with  Latin  indices,  such  as  j , running  over  the  (space)  indices 
1,  2,  3 only.  The  vector  ^ denotes  the  body  force  density,  p 


the  medium  density,  and  the  n^  are  the  components  of  tlie  normal 

* 


to  any  spatial  surface  in  the  medium.  The  fields  F^g  and  F^g 


are  given  by: 


aB 


= pVi  - "ii  ^ =1-2,3 


F, . = F. 


4j  ‘j4  P 


; j = 1,2,3 


(13d) 


F..  = p 


44 


and 


F.  . = pV.V.  - T. . ; i,j  = 1,2,3 

ij  ^ 1 J ij 


a6 


* * * 
^.1  = = Pl'j 


j = 1,2,3 


(13e) 


f = P 
44  ^ 


Similarly,  the  energy  conservation  equation,  defined  in 
of  ^ the  total  energy  density  function  for  the  medium  may  be 
expressed  as 


9x 


a 


H 


B 


where 


i.j  = 1,2,3 


aB  ^ 


- 0 ; j =.  1,2,3 


h,  -Pfv,  - Vj  ; 1 - 1,2,3 


- Pf 


and 


Hg  = (0,  0,  0,  ph) 


terms 


The  associated  condition  on  E is: 
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wi  til 


r 


K . . =0 
U 


* 

E.  . =0 
Aj 


i.j  = 1.2.3 
j = 1.2.3 


(lAe) 


E*.  = p<f ''*  - V.  T.  . + q. 
i4  1 j 1.1 


^•44  = 


i = 1,2,3 


It  is  easy  to  verify  that  the  conservation  relations  given  by  (13) 
and  (14)  are  the  same  as  those  obtained  from  equation  (12)  and 


Table  1. 
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3 . Conservation  Relations  on  : Consequences  for  Failure  Crowt h 
The  boundary  conditions  to  be  satisfied  on  E , which  may  be 
any  boundary  surface,  are: 


l»  \ < - \i'‘  "ijj  - 0 
"ij  ■» 


(15) 


(16) 


(17) 


where  the  last  equation  can  be  trivially  modified  to  include  non- 
conservative  body  force  fields  using  the  results  given  in  Table  1 
along  with  equations  (10)  or  (12) . 

These  conditions  may  be  recast  in  a more  useful  general  form 


if  one  notes  that  ||U^  "l  ]1  j;  ~ ^ uses  it  in  (10)  to  obtain 


[I  '1  "l  "l  ■ [fl  "xH, 


From  this  equation  we  obtain,  provided  F has  a nonzero  jump 
across  E : 


U.  n.  = 
1 1 


fh  "xl.-d-'i  "il, 


(18) 


A 
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The  special  condition  of  no  growth  of  the  failure  zone  is  that 

i: 


no  material  transport  occurs  across  L , in  which  case  y.  “ 0 
We  have  seen  earlier  that  this  leads  to  equation  (11).  In  this  case 
(15)  is  satisfied  identically,  while  from  (16)  and  (17)  we  obtain, 
respectively 


[tki  "ilj.  "[['kllj-  " ° 


which  are  the  well  known  boundary  conditions  at  (ordinary)  material 
boundaries.  Note  that  the  tangential  velocity  may  be  discontinuous 
on  E inasmuch  as  slip  along  E is  allowable. 

More  generally,  investigation  of  (19)  shows  that  slow  boundary 
propagation  is  characterized  by  low  values  of  the  fluxes  "ilj; 

and/or  large  values  of  EfJe  . 

We  are,  however,  most  interested  in  "fast"  processes.  From 


the  conservation  of  mass  (15),  we  have 


t 


Uf  n^ 


(19) 


In  ti’e  following  we  will  often  simply  write  ILI!  for  II]  j,  , where  it 


is  understood  that  the  jump  notation  always  applies  to  a "surface  Z 
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Let  the  unit  normal  to  Z point  away  from  the  failure  zone,  and, 
with  reference  to  figure  1,  let 

p«>  =_  P(r«>)  , p«>  . p(£«>)  . M =-  »'*’  - p"> 

where  and  denote  Z approached,  in  the  limit,  from 

region  1 and  2 respectively.  With  similar  definitions  for  the  other 
field  variables,  then  (19)  may  be  written  as 


"r  ^ ^ "I 


(20) 


where  Uj^  is  the  failure  propagation  velocity  relative  to  the 
untransformed  material.  A large  class  of  rapid  processes  of  interest 
to  us  can  be  characterized  by  small  fractional  density  jumps,  since 
we  expect  material  failure,  involving  macroscopic  and  microscopic 
fracture  and  cracking,  plastic  flow,  grain  boundary  melting,  etc  ... 
to  result  in  very  small  change  in  the  average  density.  This  kind  of 
transition  is  precisely  of  greatest  geophysical  interest  and  in  this 
case  I p|  is  very  small  and  Uj^  may  be  very  large  as  a consequence. 

We  also  observe  that  the  momentum  equation  (16)  may  be  written  in 
the  form 

- “r  fo  “i"  dh  4 - 1'’  ''J  I'*  "‘1  ■ "J 
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Eliminating  the  jump 


I"i  "il 


by  use  of  (20),  we  have 


(21) 


and  this  relation  is  exact.  When  the  density  jump  m is  small, 
then  we  have  and  so,  as  an  approximation 


Jo  \1  "r  ■ - ([' J 


(22) 


Precisely  the  same  approximate  result  is  obtained  from  (16)  if  we 

consider  U.  n.  >>  V.  n.  , which  serves  as  the  definition  of  a 
1 i 11 

fast  transition  process.  Then  n^  can  be  neglected  in  (16), 
and  we  get  (22)  directly,  with  U_  = U.  n.  . Equation  (22)  will 
be  used  extensively  in  later  sections  where  we  consider  a 
"representation  theorem"  for  the  radiated  elastic  field  associated 
with  a growing  failure  zone. 

We  note  that  typical  inferred  values  for  (shear)  stress  drops 

8 2 

and  rupture  velocities  for  earthquakes  are  of  the  order  10  dynes/cm 

5 3 

and  3 X 10  cm/sec  respectively.  With  a density  near  3.5  gm/cm  , 

then  (21)  gives  |[^kj|  cm/sec.  Thus,  (21)  implies  tangential 

particle  velocity  jumps  near  1 ra/sec  across  the  failure  boundary. 

The  energy  equation  (17)  can  also  be  put  into  a more 

appropriate  form  for  the  present  application.  In  particular, 

eliminating  [''i  "il  using  (20),  gives 


- |(v^  - q.)  n.]) 


(23) 


Here 
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1^11=  M + 1/2  [v^  vj 

where  u is  the  internal  energy  and  ((>  the  body  force  potential 
energy.  The  second  term  is  the  change  in  the  kinetic  energy  density 
of  the  material  across  the  boundary.  We  can  realistically  neglect 
the  change  in  4>  across  a failure  boundary  compared  with  the  change 
in  the  kinetic  or  internal  energies. 

The  change  in  the  internal  energy  may  be  viewed  as  the  energy 
required  for  the  transition  process  in  the  material  since  we  may 
consider  the  medium  on  either  side  of  I at  a given  time  as  two 
states  of  the  same  material.  Viewed  in  this  manner  we  observe  that, 
while  the  transition  is  irreversible,  we  can  nevertheless  consider 
u to  be  a function  of  the  ordering  in  the  material,  that  is  the 
entropy,  and  the  strain  energy  density.  In  the  transition  process 
we  would  expect  a decrease  in  the  ordering  upon  failure,  that  is  an 
increase  in  the  disorder  and  in  the  entropy,  while  the  strain  energy 
density  (i.e.,  the  recoverable  or  reversible  energy)  should  decrease 
(Here,  of  course,  the  measure  of  strain  is  relative  to  the  two 
distinctly  different  natural  states  of  the  material  in  question, 
while  the  difference  in  the  ordering  in  these  two  states  is  measured 
by  the  change  in  the  entropy.)  Therefore  we  let 

([uj  = L 


(2A) 
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with  L the  internal  energy  change  associated  witti  the  specific 
process  taking  place  in  the  material;  where  L is  used  to  emphasi;ce 
the  fact  that  the  quantity  w Is  characteristic  of  the  process 
for  the  material  and  hence  a material  property.  We  note  that, 
because  of  the  definition  of  ([uj]  = u^^^  - u^^^  , we  expect  the 

change  in  u across  Z due  to  disordering  of  the  material  to  be 
negative,  while  the  part  of  the  change  due  to  the  strain  energy 
positive.  Thus  the  total  change,  being  the  sum  of  these,  may  in 
principle  be  positive  or  negative.  We  will  show,  however,  that  for 
a spontaneous  process  to  occur  along  with  rapid  growth  of  the 
transition  zone,  then  it  is  generally  necessary  that  L be  positive. 

In  particular,  if  the  spatial  heat  flux  jump  across  Z is  neglected 
(i.e.,  if  the  change  in  the  thermal  gradient  and  possible  changes  in 
the  thermal  coefficient  across  Z are  neglected  as  small  compared 
with  the  other  effects),  then  L must  always  be  positive.  In 
general,  for  failure  processes  we  expect  terms  like  tj^j  in 

(23)  to  be  much  larger  in  magnitude  than  the  term  ||qj^  nj  , and 
this  size  ordering  also  leads  to  the  requirement  that  L be  positive 
for  non-zero  growth  rate  of  a spontaneous  failure  zone.  A requirement 
that  L be  positive  simply  means  that  the  change  in  the  strain  energy 
must  generally  be  larger  in  magnitude  than  that  part  of  the  internal  energy 
change  due  solely  to  disordering.  This  requires  relatively  high 
(nonhydrostatlc)  strain  energy  for  such  a process  to  occur. 

To  show  that  these  statements  follow  from  the  condition  (23), 
we  can  rewrite  it  as 
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K 'k 

“k  J!l.kj  .. 

P^L  + 1/2 

I'k  vl) 

where  we  consider  the  case  In  which  the  denominator  is  non-zero, 


of  course.  Here  we  have  suppressed  the  superscript  on  the 
,(1) 


density,  p , and  simply  written  p to  denote  this  quantity. 
Using  (21)  in  this  equation,  after  expanding  both  I''k  'kj  ■"'* 

[wl  • 


we  have  the  results: 


Ur=  - 


h 

1 ■ 

-I 

Ik  \ 

hk  - k \ 

P I 

. + 1/2 

M'-kl 

p 

Using  (21)  again,  in  either  of  the  expressions  for  Uj^  in  (26), 
to  eliminate  terms  in  Set* 


n2  _ n _ tjijJ  = 


R pL  R 


2 p^L 


The  roots  of  this  quadratic  in  Uj^  are  real  provided 


ti2 


+ 2L 


‘'k 


> 0 


(25) 


(26) 


(27) 


(28) 


When  this  condition  is  not  met,  then  Uj^  has  complex  (or  imaginary) 
roots  and  this  means  that  initiation  and  growth  of  the  transition 
process  cannot  occur. 


no 


Wlien  the  st 


then  ||t 
are  all 


k ‘^k 


> 0 

positive . 


resses,  and  hence  the  tractions,  drop  upon  failure, 
, and  the  scalar  terms  involving  the  jumps  on  T. 
In  this  case  we  must  have 


(28a) 


As  was  already  noted,  for  a failure  process  we  expect  the  heat  flux 
term  to  be  much  smaller  than  the  traction  term  so  that,  for  all 
practical  purposes,  we  must  have 

L > 0 


If  we  neglect  the  heat  flux  term  completely  on  the  grounds  of  its 
relative  size  for  a rapid  spontaneous  failure  process,  then  (27) 
gives 


(29) 


Alternately,  under  these  circumstances 


(30) 


We  note  that  rapid  growth  of  a tectonic  failure  zone  Involves 

Uj^  = 3 X 10^  cm/sec  and  traction  jumps  at  the  front  of  the  failure 

8 9 2 

boundary  which  appear  to  be  of  the  order  of  10  to  10  dynes/cm 


Ill 


(e.g.  , Archambenii,  1976).  This  implies  that  the  internal 

4 6 

energy  change  I.  is  in  the  range  from  10  to  10  ergs/gm. 

We  can  also  express  L in  terms  of  the  particle  velocity 

and  traction  changes  across  X when  we  neglect  the  heat  flux  term. 

2 

That  is,  eliminating  from  (30),  using  (21),  gives: 


L - in  v,^ 


'1 


The  ratio  of  the  quadratics  in  traction  in  (31)  will  be  near  unity 
when  the  hydrostatic  stress  is  maintained  after  failure  (i.e.,  small 
specific  volume  change  and  density  change)  and  when  the  deviatoric 
stresses  nearly  vanish  upon  failure — that  is  when  there  is  nearly  a 
complete  loss  of  shear  strength  upon  failure.  It  is  quite  possible 
that  this  situation  would  prevail  for  many,  if  not  most,  earthquakes. 
If  this  is  the  case  then 


L ~ 1/2 


I"  kD  [''kj 


(31a) 


With  L in  the  range  10^  to  10^  ergs/gm  then  this  implies  a magnitude 
of  from  1 to  10  m/sec  for  particle  velocity  jumps  across  the  failure 
surface  and  conversely. 
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If  the  heat  flux  term  can  be  ne{>.lected,  that  is  if 

I d^k  "kj  ^ ^ [[.'^k  '^kj  relation  (23)  or  (25) 

may  be  expressed  as 


If  we  use  the  expression  (31)  for  L in  this,  then  we  find  that  this 
expression  of  energy  conservation  can  be  put  in  the  form: 

In  view  of  the  relation  (21)  for  , we  see  that  this  equation 

is  satisfied  when  momentum  is  conserved. 

conservation  of  energy  is  insured  when  equation  (31)  (or  (31a))  is 
satisfied  along  with  momentum  conservation — which  is  expressed 
by  equation  (21).  Consequently,  the  equations  for  failure  zone 
growth  seem  to  be  well  approximated  by  the  relations 


k "k 


J 


then 


Therefore  this  shows  that  if 


I (32) 
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i 


; 


If  we  treat  L as  a material  property,  then  the  second  of  these 
relations  can  be  viewed  as  a condition  for  failure  grouth--that  is 
when  the  particle  velocities  and  tractions  across  the  failure 
boundary  are  such  that  this  relation  is  satisfied  for  a value  of 
L appropriate  to  a specific  failure  mode  of  the  material,  then 
growth  may  occur;  with  the  rate  of  growth  then  given  by  the  first 
relation  in  (32).  In  at  least  some  instances  of  rapid  failure 
the  shear  tractions  would  be  very  small  within  the  failure  zone, 
and  then,  as  previously  noted,  the  second  relation  reduces  to 

Thus  we  may  view  the  second  relation  in  (32)  as  a dynamic 
failure  condition  and  the  first  relation  as  an  equation  for  U_  , 
the  dynamic  growth  rate. 

The  dependence  of  the  rupture  growth  rate — or  rupture 
velocity  — on  the  rheological  properties  of  the  material  can  he 

inferred  from  the  first  relation  in  (32)  if  we  specify  appropriate 
constitutive  relations  for  the  material  in  its  two  states,  before 
and  after  failure  has  taken  place.  A constitutive  relation  of 
sufficient  generality  to  describe  the  local  behavior  of  the 
material  after  failure  is  of  the  form  (e.g. , Fung  1965) 


T.  . 
ij 


■E/.  “ 


(t-T) 


9e 

r p + r * 
i.ikP.  k^  ^ijkf  8t 


kf. 


(33) 


I 

I 


I 
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where  the  term  involving,  the  sum  over  Q arises  from  relaxation 

processes  within  tlie  material,  of  which  there  may  he  N types  each 

denoted  by  an  index  a . Here  is  a relaxation  frequency 

(reciprocal  relaxation  time),  and  the  object  is  a time 

ijkil 

independent  relaxation  modulus  (or  "relaxation  strength")  tensor. 

Relaxation  terms  can  be  shown  to  result  from  a variety  of 
processes,  including  the  movement  of  interstitial  atoms,  and  vacancies, 
twining,  chemical  reactions,  crystalline  thermal  currents  and 
dislocation  motion.  Descriptions  of  typical  physical  processes 
of  this  type  for  solids  and  liquids  are  given,  for  example,  by 
De  Groot  and  Mazur  (1969)  and  Mason  (1966). 

The  second  term  in  (33)  is  a linear  elastic  term,  while  the  last  term 
is  a simple  linear  viscous  term. 

The  strain  is  defined  as  an  infinitesimal  strain  measured 

from  the  completely  relaxed  state^  of  the  material  aftev  transition. 

Thus,  in  the  infinitesimal  strain  approximation  which  we  will  adopt 
for  the  material  after  transition: 


. (2)  . (2)) 

= <2)  - w,  I 

®kJ!.  “(k,Ji)  I ax.  ax, 

it  K 


(2) 

where  u,  is  a displacement  from  the  relaxed  reference  state, 
k 


(34) 


t 


By  relaxed  state  we  mean  the  equilibrium  configuration  taken  by  a 
piece  of  the  material  removed  from  the  surrounding  medium,  so  that 
it  is  free  from  external  forces  and  surface  tractions. 
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A similar  constitutive  relation  liolih:  for  the  material  outside 
the  failure  zone  or  for  the  material  before  failure.  lloweve. , the 
relaxed  state  of  this  material  is  different  from  that  tor  the  failed 
material,  so  that  the  infinitesimal  strain  is  defined  as 


= (1)  =1/91 

“(k.t)  “ 3x^  9x 


with  Uj^  the  displacement  measured  relative  to  the  relaxed  state 
of  the  unfailed  material. 

Further,  while  the  form  of  the  constitutive  relation  may  be 

(a) 

the  same  for  the  material  in  the  two  states,  the  moduli  D , „ 

ijk£ 

^ijkil  ^ijk£  very  different.  For  the  unfailed 

material  we  can,  in  the  present  context  at  least,  neglect  the 
(ct) 

anelastic  moduli  D. „ and  viscous  term  C! „ relative  to  the 
ijkJi  ijkJt 

linear  elastic  term  C...^  for  the  unfailed  material.  Thus  for 

ijk£ 

the  region  outside  the  failure  zone  we  may  use  a constitutive 
relation  of  the  form 


T = r e 

ij  iik£  k£ 


with  defined  by  (35).  In  any  case,  however,  (33)  can  be  used 

to  represent  the  form  of  the  constitutive  equation  in  both  the 
failed  and  unfailed  material,  with  the  coefficients  and  strains 
taking  on  different  values  in  the  two  zones. 
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We  can  write  (33)  in  the  compact  operational  form 


where  Integral  operator 


“..u  = 5 J- 


t -!!„<t-l:') 


+ «2g 


with  denoting  the  Kronecker  delta,  while  6(t-t')  is  a Dirac  delta 

(6) 

function  and  6,(t-t')  its  derivative.  Here  m, . is  a time  independent 
1 i ] kJ- 

modulus  tensor,  equal  to  the  elastic  and  viscous  tensors  for  3 = lf2 
and  equal  to  the  set  of  relaxation  moduli  tensors  for  3 = 3,4  ... 

We  can  also  write  (38)  as 


^'ijkSi  " S "’JjkSi  ^3 


where  I_  is  the  time  dependent  integral  operator  in  (38). 
3 

Now  the  momentum  condition  in  (32)  can  be  written  as 


■ ■ M ■ ■ 'e J 
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Without  undue  loss  of  generality,  we  will  consider  the  isotropic 
case.  This  gives 


ijk«. 


= A-  6.. 
S i: 


Ik 


^jk 


<5.,) 

ik  j£. 


for  the  modulus  tensor. 

If  we  now  define  a signal  velocity  Cg  to  be  associated  with 
the  rate  at  which  changes  in  the  displacement  fields  are  propagated, 
in  the  sense  that  displacements  are  causal  with  respect  to  such  a 
velocity,  then 

(41) 


where  n = 1,2  for  the  two  material  states,  and  T 
Then  we  have: 


t - r/C^  with 


r E (x.  x.)l/2 
J J 


8u 


(n) 


dx„ 


JL 

c. 


(?) 


8u 


(n) 


3t 


Since  we  can  use  the 

infinitesimal  strain  approximation  for  the  material  in  its  initial 
(unfailed)  and  final  (failed)  states,  then  the  reduced  time 
derivative  of  u^  is  also  the  particle  velocity.  Thus 


y(n) 


3x 


I 


(A2) 
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Using  (42)  in  tlie  mointintum  condition  yields 


P "r  li\l  - s 

P 


(43) 


for  the  isotropic  case.  The  magnitude  of  U„  is  thus  directly 

proportional  to  the  changes  in  the  material  moduli  upon  failure. 

For  example,  suppose  the  material  has  the  properties  of  a perfect 

fluid  after  transition.  If  the  jump  in  the  normal  component  of  velocity 

across  the  transition  boundary  is  small  compared  with  the  jump  in  the 

tangential  components  of  the  velocity  (i.e.,  small  changes  in  density 

and  compressibility)  then  we  have  as  an  approximation  to  the  maximum 

magnitude  of  U : 

R 


considering  the  failure  process  to  be  driven  by  a shear  wave,  so  that 


Cg  = Vg  with  Vg  = /y/p  the  shear  velocity  in  the  medium  before 
failure  and  since  y = 0 for  an  ideal  fluid,  then: 


V 


S 


In  conclusion  then,  it  is  evident  that  the  form  of  (43)  indicates 
a dependence  of  the  failure  rate  on  the  transition  process  in  terms  of 
the  changes  in  material  moduli.  In  general  we  see  that  consideration  of 
an  expression  like  (40)  is  required  in  order  to  obtain  a precise  estimate 
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of  the  magnitude  and  space-time  dependence  of  U„  . Nevertheless, 
simplified  expressions  such  as  (43)  and  approximate  results,  such 
as  that  for  the.  ideal  fluid  transition,  are  quite  useful  and  can  provide 
a basis  for  interpretation  of  observed  failure  rates. 


120 


4 . Linear Lzat ion  of  tlie  C oupled  Conservation  Egua ions 

We  will  now  focus  our  attention  on  the  medium  outside  the  failure 
zone  and  express  the  conservation  relations  in  this  region  in 
linearized  form,  using  the  usual  approximation  of  infinitesimal 
strain  for  this  region.  This  will  give  linear  equations  of  motion 
for  this  part  of  the  medium  when  we  take  tiie  material  to  be  elastic 
and  specify  the  rheology  according  to  (34). 

We  note  that  we  can  also  describe  the  dynamical  behavior  of 
the  medium  within  the  failure  zone  itself  after"  failure  using  linear 
theory,  if  we  adopt  (33)  as  an  adequate  description  of  the  rheology. 

In  this  case  we  would  simply  repeat  the  arguments  leading  to 
linearized  equations  and  obtain  an  integral  Green's  function 
solution  in  a manner  analogous  to  tliat  for  the  region  outside 
the  failure  zone.  The  two  integral  representations  of  the  medium 
response  would  then  be  connected  by  the  boundary  conditions  on 
the  failure  surface — as  given  by  the  equations  (32)  in  tlie  previous 
section.  In  this  development,  however,  we  will  only  consider  the 
exterior  region.  This  does  not  result  in  any  loss  of  generality 
in  the  analytical  description  of  the  radiation  field  since  the 
dynamical  beiiavior  of  tlie  in..erior  medium  manifests  itself  directly 
in  tlie  boundary  conditions  for  the  exterior  problem,  and  these  can  be 
treated  in  a formal  manner  with  the  tractions  and  tiie  particle 
velocity  of  the  interior  zone  material  at  the  failure  boundary 
being  left  completely  arbitrary.  Tlie  resulting  general  integral 
rcpresentaf ion  for  the  exterior  radiation  field  can  then  be 
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i 


t 


i 


evaluated  for  particular  cases  when  the  rheological  properties 
of  the  material  after  transition  are  independently  specified  . 

Then  the  resulting  dynamical  variation  of  tlie  boundary  surface 
tractions  and  the  particle  velocity  are  determined  by  either  solving 
for  the  interior  deformation  fields  using  the  Green's  function 
representation  or  by  choosing  a simple,  yet  adequate,  rheological 
description  for  which  the  dynamical  behavior  in  the  interior  is 
known  with  sufficient  accuracy  a priori.  An  example  of  the 
latter  would  be  the  assumption  of  melting  and  the  production  of 
an  essentially  inviscid  interior  fluid,  so  that  the  shear  tractions 
in  the  interior  would  vanish  and  the  particle  velocity  would  only 
reflect  any  accompanying  density  change. 

To  treat  the  exterior  region  we  first  observe  that,  in  the 
present  context,  it  may  be  viewed  as  containing  an  elastic 
material.  In  this  case  the  energy  equations  are  trivially 
satisfied  in  the  medium  away  from  the  failure  boundary  if  we  neglect 
heat  transport  and  transport  terms  in  the  energy  equation.  On  the 
failure  boundary,  however,  energy  is  absorbed  by  the  failure 
process  and  energy  conservation  is  expressed,  approximately,  by 
the  second  relation  in  (32). 

The  linearized  equations  of  motion  for  the  exterior  medium  can  be 
written  in  the  form 


9 

9t 


1 3 u.\ 

9 u 

( n — ^ = -5_ 

f Ji 

8t  / 9x. 

J 

iikP.  9xj^ 

+ p b . 
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(4  4) 
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whert-  IK  is  the  elastic  displacement,  and  both  the  density  p 

and  the  elastic  tensor  C..,„  are  treated  as  functions  of  the  spatial 

1 j kx. 

coordinates  alonet  In  obtaining  (43)  we  have  used  the  symmetry  relations 
(e.g.,  Landau  and  Lifschitz,  1965): 


*^ijkZ.  ikX,  *^ij£k  ^k£ij 


Because  of  the  form  of  (44),  we  can  write  these  linearized 
equations  of  motion  in  the  same  four  dimensional  form  as  was  used 
to  express  all  the  conservation  relations  earlier.  Specifically, 
we  define  an  "elastic-inertial  tensor"  as: 


^ijk£ 


ajB.Y.S  = 1,2,3 


^aBY6  • i '^i4k4  ^4i4k  ^ “^ik  ’ 


C „ = 0 ; otherwise 

aSyo 


where  C . . is  seen  to  have  the  symmetry  properties 
apyo 


*^a66Y 


tin  what  follows,  we  use  the  common  notation 


k,£  3x. 


for  partial  derivatives. 


u = — > etc. 

k,£m  9x„  3x 
Jt  m 
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with  tl\e  Creek  indices  ranpjng  from  1 to  A.  If  We  also  define  tile 
spacclike  fields 


u - (u  , u , u , 0) 
a 12  3 


f H (b. , b„,  b„,  0) 
a 12  3 


(A7) 


then  (A4)  can  be  written  in  the  form 

V.,6  = 


a6 


(48) 


Alternately,  in  a form  convenient  for  a solution  for  the  displacement 
field  , the  equations  of  notion  may  be  expressed  in  the  operator 
form: 


L u = pf  (49) 

ay  Y a 


where  tlie  "elastic  operator"  is  given  by 


L 

ay 


(50) 


4* 

Kupradze,  1963,  defines  the  "elastic  operator"  in  an  analogous 
fashion  in  the  frequency  domain. 
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The  boundary  conditions  can  also  lie  expressed  in  terms  of 
the  elastic-inertial  tensor  defined  in  (‘^*1)  and  the  spacelike 
displacement  field  u.^  if  we  define  a four  dimensional  space-time 
"normal"  as 


= (Oj^,  n^,  n^,  - U 'n) 


(51) 


with 


where  n.  is  the  regular  spacelike  normal  to  any  spatial  surface 

■k  ^ 

in  the  medium,  in  particular  the  surface  E , and  U ’n  = U is 

— K 

the  projection  of  the  spacelike  relative  velocity  vector  for  the 
surface  E along  the  normal  to  the  surface.  Now  we  observe  that 
in  view  of  the  properties  of  the  elastic-inertial  tensor 


aByS  '^y,6 


kilmn 


u 

Tn,n 


P 


k,4 


* 


1 


or 


k 


- 


but  the  right  side  expresses  conservation  of  momentum  on  E , 
is  verified  from  (16),  and  vanislies.  Thus  we  have  that 


as 


^'3  Jl,  = « 


for  all  a and  3 . 

Thus  the  linearized  equations  for  tlie  exterior  region  are  most 
naturally  expressed  in  a four  vector  form,  as: 


L u = T„„  = pf 

ay  y a3,3  a 


K3  "3  1 = 


with  the  inertial-stress  tensor  given  by  (48)  and  where  the 

elastic  operator  is  defined  by  (50)-  In  the  following  we 

will  linearize  the  boundary  condition  in  (52)  to  the  extent  that 
we  will  use  instead  of  , which  amounts  to  the  neglect 

of  Vy  compared  with  , as  is  justified  for  failure  processes. 
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5 . Green's  Tensor  Equations  for  the  Radiation  Flolcl  _i n_ 

Elastic  Zone: Elastodynanlc  Representation  Tj>e:£.tema 

The  different i>al  operator  L ajiplies  to  space-1  U'.e 

ay 

tensors  {i.e.,  tensors  without  time-like  components)  with  the 
property , 


=”a6...X 


The  contraction  L w . , is  another  space-like  tensor. 

ua  aS.  . . X 

The  operator  has  been  defined  in  a four 

dimensional  space-time,  and  it  is  useful  to  consider  an  arbitrary 
four  volume  in  this  space  over  which  an  inner  product  of  space- 

like tensors  is  defined.  Thus,  given  two  space-like  tensors 

e and  h , we  define  an  inner  product  in  il  by 

^a...X  6...y 


t 


Here  (g  , g . ) is  positive  and  vanishes  only  when 

g , H 0 . Let  C*^(x,  x')  denote  a tv/o  point  tensor,  a function 

a.  . . A ry  ’ 

of  the  independent  coordinates  x and  2i'  in  51  , then  ^^^e  define 

the  inner  products 


(G 


(G 


6 

a 


u )-'  - f G^(  X,  X')  u (X') 

a ^ J a - - 

n 

U )-  r.  /'g^(  X,  X')  u ( x)  d^x 
a f,  y a n - 

Q 


We  will  use,  as  a notational  convenience,  both  superscript  and 
subscript  indices  for  the  tensors  appearing  here.  However,  all 
tensors  are  Gartesian  and  there  is  no  distinction  lietwee>n 
covariant  or  contravar iant  tenfiors. 
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which  give  a vector  funecion  of  x in  tlie  first  ease  or  a vector 
function  of  x'  in  the  second.  If  C is  symmetric  in  x and 
x'  , then  the  two  inner  product  forms  are  entirely  equivalent. 
The  Green's  functions  associated  with  the  operator  L are 

ay 

two  point  tensors  tliat  are  of  special  importance  in  the  present 
development.  In  particular,  we  define  the  two  point  Green's 
tensor  associated  with  to  be  given  by 


L G®(  X,  X ) = X,  X ) 

ay  y 0 a ~ —Q 


where  the  space-like  tensor  is  a fundamental  solution  of  the 

operator  with  a pole  at  ; with  the  space-like  tensor 

defined  to  be  the  generalized  delta  function: 
a 


Here  6(  X“  2^^)  is  the  (four-dimensional)  Dirac  delta  distribution. 

Clearly  G defined  by  (5J)  will  have  a singularity  at  x = x due  to 
Y ~ ~o 

the  presence  of  the  Dirac  delta  distribution  and  therefore  will  satisfy 

(5J)  in  a distributional  sense  (see  for  example:  Stakgold,  1968).  In 

the  present  context  G^(  Jc,  x^)  can  be  viewed  as  the  y component  of 

the  displacement  field  in  the  continuum  at  21  due  to  an  impulsive 

force  in  the  B direction  located  at  x . Consequently  x can  be 

~o  o 

thought  of  as  the  source  point  and  21  as  the  receiver  or  observ'er  point. 

With  these  definitions  of  the  inner  product  in  D for 
space-like  tensors  we  can  generate  integral  equations  that 


128 


are  ec|iiivalouL  to  tlu;  differentia]  iMiuationn  and  anatii' iat  C'd  bound  iry 
conditions  w’nic.h  specify  the  displacement  field  in  a continuum. 

To  do  so  we  will  make  use  of  the  special  propi-rties  of  the  Green's 
function. 

In  particular,  the  Green's  theorem  in  12  is  obtained  by  consideriiig 
the  inner  product 


(1-  w.  v) 


12 


(54) 


where  the  operator  L is  simply  the  matrix  operator  with  components  given 
in  (49),  while  w and  v are  fields  satisfying  eitlier  of  the  operator 
equations  (48)  or  (53)  and  hence  may  be  single  or  two  point  space-like 
tensors  of  first  or  second  order. 

Now  using  the  representation  of  the  components  of  L given  in  (49) 
we  have 


(Lw 


/ 


d^ 


(55) 


A result  similar  to  the  ordinary  Green's  theorem  is  obtained  by 
first  observing  that  die  integrand  can  be  expressed  in  terms 
of  an  identity  as: 


V 

a 


(v  C „ . w ) 
a a3Y'5  Y."  ,3 


-V  „C„^w 

a, 6 aSYiS  Y»o 


However  the  final  term  can  be  rewritten  by  using  the  symmetry 

properties  of  C , so  that 

' apiyo 


I 

i 


I 

i 

i 

i 

1 


I -• 


■V  „C„,W  c=w(C  v t)D-(w  C „ , w 
u,B  aBY‘5  Y>5  a oBy^  Y>S  ,3  u ctBYiS  Y»o  >B 


Using  these  results  in  (55)  we  get 


adY^  Y.iS  .3 


■*■  / [(v  C w r-w  C .v  .)  ]dx  (56) 

/ a aSYS  y,&  o uBy^  YjS  >3  ^ ' 


We  can  now  define  a formal  adjoint  operator  L , where 


/ = (C  — — ) 

aY  3Xg  aBYiS 


are  the  matrix  elements  of  L and  further  define  the  bilinear  concomitant  J. 

t 

of  w and  v where 


Jn  = V C„oW  ^-W  CorV  ^ 

3 a a3Y<5  Y>o  « a3Yo  y,o 


in  order  to  write  (56)  in  the  form 


(Lw  .V  )|  = ( w,Lv  )^  + f 


This  is  the  generalised  Green's  theorem  for  the  operator  L in 
We  observe  from  comparison  of  (57)  with  (49)  that 
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L 


♦ 

L - L 


and  hence  that  L is  formally  self  adjoint. 

In  deriving  (59)  we  have  considered  the  operation  of  L on  fields 
specified  at  x and  h.ave  formed  inner  products  on  Q with  respect  to 
the  coordinates  x . We  may  form  inner  products  with 

respect  to  fields  specified  at  the  coordinates  , which  in  the  case 

of  two  point  tensor  fields  such  as  ( 2<.  2^)  means  integrating  over 

the  "source  coordinates"  x rather  then  the  "observer  coordinates"  x . 

— o — 

In  formal  terras  we  may  generate  a completely  parallel  result  to  (59)  . 
Consider  the  operator  defined  by 


(60) 


in  component  form.  Now 


L u 

ay  i 


3x, 


a3y<5 


Su 

3x 


J 

■^6  ^ 


is  just  the  operation  of  L on  u at  x . Here  of  course  u is 

- -o  Y 

expressed  as  a function  of  the  coordinates  x . In  case  the  operand 


is  a two  point  tensor  G 


3 


G®  ( X ; x) 

ay  Y — o — 


2io’  * 

o 

dx. 


• 4 


F 
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Clearly  we  can  nuw  form  the  inner  product 


//(O) 

(i-  y.. 


V L w d X 

a ay  y 


By  exactly  the  same  formal  manipulations  used 

to  obtain  the  previous  result,  we  have  the  complementary  Green's  theorem 
in  fl  for 


(o) 

^ ‘'r-T  tr  \ — 


, z) 


where  is  the  adjoint  of  and  where 


l_(o)*  ^ ^(o) 


so  that  is  formally  self  adjoint,  as  is  obvious  since  L was 

self  adjoint.  Here  also  is  formally  identical  with  , with 

the  coordinates  x replaced  by  the  coordinates  x° 
p p 

If  w and  V are  one  point  tensor  fields  then  (59)  and  (62)  are 
clearly  identical  results.  However  if  one  or  both  of  w and  v are  two 
point  tensor  fields  of  x and  21^  » then  (59)  yields  relations  between 
function  of  x^  while  (62)  yields  relations  involving  functions  of  x 
which  may  be  different. 

The  integral  relations  obtained  take  a physical  meaning  when  the 
fields  Z and  Z are  specified.  In  particular  we  will  take  w to  be  the 
displacement  field  u satisfying  the  equations  of  motion  and  boundary 
conditions  for  the  region  fi  surrounding  a closed  failure  region.  In  this 
region  we  take  the  failure  zone  to  have  a volume  with  a closed 

boundary  surface  T,  while  tlie  rcmiai.ning  spatial  volume  is  denoted  by 
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with  an  external  boundary  S , as  shown  in  Figure  2.  Thus  tlie 
boundary  of  is  Y,  ® S while  the  boundary  of  is  Y.  . The 

space-time  volume  fl  is  the  four  volume  occupied  by  the  medium  externa] 
to  E as  shown  in  Figure  2.  The  external  boundary  of  fl  is  generated 
by  parallel  displacement  of  S in  the  direction  of  the  (time)  axis, 

since  S does  not  change  with  time.  The  internal  boundary  of  Q , however,  has 
generators  that  are  not  parallel  to  the  x^  axis,  reflecting  the  growth 
of  the  failure  volume  with  time. 

In  addition  to  the  space-like  boundary  conditions  expressed  by 
the  boundary  conditions  in  (52),  which  apply  on  the  spatial  boundary 
Y , there  may  be  conditions  on  the  fields  that  apply  on  the  time-like 
boundary  of  fl  . These  may  be  connected  with  the  space-like  conditions 
in  (52)  in  that  they  may  result  from  the  change  in  E as  a function  of 
time,  and  in  any  case  the  conditions  along  the  time-like  surface  of  (2 
would  correspond  to  "generalized  initial  conditions".  These  general 
time-like  conditions  are  in  fact  already  contained  in  the  boundary 
conditions  (52)  provided  we  recognize  that  E is  a function  of  time 
and  that  the  fields  may  be  discontinuous  in  time  as  well  as  in  space. 

We  will  later  show  how  the  time-like  discontinuous  behavior  gives  rise 
to  relaxation  phenomena  that  account  for  the  stress  wave  radiation 
accompanying  the  growth  of  the  failure  zone  in  an  initially  stressed 
medium. 

To  complete  the  description  of  the  physical  problem  we  consider  the 
impulse  response  of  the  medium  to  a simple  delta  "function"  point 
source.  Tt’.is  response  is  of  course  provided  by  the  Green's  function 


appearing  in  equation  (53),  the  role  of  Green's  integral  theorem  being 
simply  to  provide  the  means  of  superposing  or  adding  together  weighted 


Figure  2:  Geometry  of  the  four  volume  Q with  x.  and  spatial  coordinates, 

X the  time  coordinate.  A spatial  section  through  the  four  dimensional  volume 

at  time  t,  is  the  ordinary  spatial  volume  V, (t,),  external  to  the  failure  region 


Cii'fu’s  functions  in  such  a way  as  to  reprosunt  the  coniplox  si'uri-u  wi- 
are  cicalin;’  witli. 


r 


I 


i:i4 


Therefore  we  will  Lake  the  field  V appearing  in  the  Green's 
theorem  formulae  (59)  and  (62)  to  be  the  two  point  tensor  field  x^) 

We  shall  require  that  causality  is  satisifed,  that  is,  (e.g..  Horse  and 
Feshbach.  1953) 


(X 


*o)  = 


for 


(63) 


and  hence  that 


^3  / . o o.  _ . o 

^4^  ^Y  ^k 


(64) 


In  addition  to  the  differential  operator  equation  and  the  causality 

D 

condition,  is  only  fully  determined  by  specification  of  boundary 

a 

conditions  to  be  fulfilled  and  the  choice  of  boundary  conditions  appropriate 
for  G*^  depends  upon  the  physical  problem  to  be  solved.  In  the  present 
case,  the  problem  is  represented  by  the  operator  equations: 


Lu  = pf  ; xef2 

B U = b ; Xc3P. 


(65) 


The  second  equation  simply  represents  the  boundary 
in  compact  operational  form.  In  particular  B is 
components 


conditions  in  (52) 
a matrix  operator  with 


Ruby'S 


3x . 
0 


I 

i 


I 
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so  that  5 Op  • Further,  on  Z , is  the  four  vector  equal  to 

the  value  of  i - ii„  arising  from  the  action  of  the  material  within 
an  n 

the  failure  zone  enclosed  by  Y.  , which  we  may  take  as  independently 
specified.  For  a physical  system  governed  by  operator  equations  of 

g 

this  form,  then  G S j taken  to  satisfy  the  homogeneous 

system  (e.g.,  Stakgold,  1968);  | 


LG 


= A ; Xefi 

= 0 ; Xean 


(66) 


I 


with  L and  8 being  identical  operators  in  (65)  and  (66)  . We  note  that 
the  Green's  theorems  expressed  by  (59)  and  (62)  involve  the  adjoint 
operators  L*  and  ^ Consequently  we  also  define  a system 

adjoint  to  (66)  as 

L*  G*  = A ; xeP. 

8*  G*  =0  ; xeSP 


! 

I 

I 


where  G is  defined  as  the  adjoint  Green's  function,  L*  the  operator 
adjoint  to  L and  8*  Q*  = 0 as  adjoint  boundary  conditions  defined  to 
be  such  that  ; 


ac.G-)”-  <G.  CG-)“ 


(08) 


I. 


I 


I 
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ll 


Rfferriiig  to  the  results  (‘ly)  - (g21  w*?  observe,  for  the  system  represeateO 
by  the  clLfferciitLal  operator  and  boundary  conditions  of  (65)  and  the  associated 
systems  (66)  and  (67),  that  8*  is  determined  by  the  condition  that  the 
bilinear  concomitant  of  G and  G vanish. 

Noting  that 


‘■ay  ''■y  < > • V ’ ^ 


•■ayS  ■ ''a'*^  «l) 


with  and  arbitrary  source  points  in  Si  , and  inserting  these 

expressions  in  the  integral  relationship  (68  ) involving  G and  G*  , 
we  iiave,  using  the  properties  of  the  delta  functions  : 


G,^"(  X_;  X ) = ( X ; Xy)  ^ 


Y o 


ft*  ft 

g'^  (X  ; X„)  = G*"  ( X^;x  ) 
Y o Y o 


(6^11 


which  is  the  reciprocity  relation.  Thus,  from  (64)  and  (69)  it  follows 


that  the  causal  Green’s  tensor  iias  the  adjoint  given  by 


o o 


-I;  -x^  x,.,  -xj 


Y " k’  4’  k’  Y k’  A’  ""k’  A' 


(70) 


3* 


Gy  (Xj^,  x^;  x^^,  x^^)  = G (Xj^,  -x^;  x^,  -x^)- 


(o) 


Now  consider  the  operator  I,  ' , defined  to  act  on  the  source 


coordinates  X . Clearly,  i\^torchanging  the  roles  of  X and  X^  (66)  gives 


'S’  ■ “!  "‘o'*'  “ S 


B ( X ; X ) = 0 

ay  Y o 


By  the  reciprocity  relation  (69),  then 


l(o)  X ; X ) 

OtY  Y ’ o 


A^(  x;  X ) 
a o 


Similarly, 


B (X  ; X ) =0 

ay  Y o 


, (o)*  pB  , . 

ay  Y ^ ’ *o 


6„  G«  ( X ; X ) - 0 

ay  Y o 


Now  returning  to  the  Green's  theorem  expressed  by  (62),  where 
the  integration  is  over  the  source  coordinates  X^  , and  using  G for 
V and  U for  w , we  have 


(1^°^  u,G)*°  = (u,  1^°^'"G)* 


ri 


,G)  d^x°  (73) 


where  is  seen  from  (68)  to  be 

p 


(X  ; x^)  T .(  X ) - u ( X^)  G^'  ( X ; X ) 

6a  oapo  aoa3  o 


(74) 
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Now  formal  application  of  Clio  divcryonco  Lliiorom  to  the  second 
integral  on  the  right  gives: 


- u G^'J 

d\«  = [ 

f’l*-  pb 

1^0 

'r. 

o 

a rtB 

a ad 

J 

dil 

a a3  a nS 

.3  o 
d X 


(77) 


The  integral  over  the  boundary  9::^  however  takes  a special  form 

because  some  of  the  fields,  namely  u and  G , are  purely  space— 

Cl  a 

like,  while  both  t „ and  have,  by  construction,  tine-like 

ap  ap 

components.  Hence,  we  can  give  more  explicit  form  to  the  "surface" 

integral  in  (77)  by  applying  the  divergence  theorem  to  the  expanded  form 

of  the  integrand,  where  and  are  written  out  in  terms  of  the 

purely  space-like  fields  used  in  their  construction.  From  the  definitions 

of  T and  G'\,  we  have 
ap  a£i 


f 

9 

G T ^ - u 0 

1 

n 

9x° 

ot  a3  a otS 

B 

A o 

d x 


9u . 


J 


3u . 
1 


3G 


kHij  Bi 


9x 


ac"  n 

s 


A o 

d X 


where  the  latin  indices  take  only  the  values  1,  2,  3.  Regrouping  the 
terms  and  noting  that  the  ordinary  elastic  stress  is 
similarly  that  clastic  stress  associated  with  the 

displacement  C™  , we  have: 
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f _J_  , g"  1 = f 

I o \ a a ap  J J 

dX  L o 


0 '"b 


0 


,,m  r 


] d'^o 


n ““'A 


/ ^^1 
r™  f „ 

d-  H 


.4  o 
d X 


(78) 


Now  we  observe  that  the  first  integral  involves  a purely  spatial 
divergence  of  a space-like  field,  while  the  second  integral  involves 
the  purely  time-like  part  of  the  (four  dimensional)  divergence  of  a 
spatial  field.  Thus  we  integrate  over  only  spatial  coordinates  in  the 
first  and  over  the  time  coordinate  in  the  second,  in  the  latter  case 
observing  that  fi  is  an  explicit  function  of  time  since  the  surface 
T.  is  expanding.  In  particular,  taking  the  initiation  of  the  process 
resulting  in  a dynamic  field  to  begin  at  ^4  = = 0 and  noting  that 

g"*  is  causal,  so  that  C™  = 0 when  , then 

k K 


kX. 


“A 


A o 

d X 


/ [ 


„m 


- uj^  G 


u]", 


da 


( 79) 


by  an  ordinary  application  of  the  divergence  theorem.  Here  is 

the  boundary  of  the  spatial  volume  at  time  t , and  n^  is  the 

normal  to  this  spatial  boundary,  as  siiown  in  Figure  2.  Also,  t denotes 

+ 

t + e , with  f.  > 0 and  infinitesimal.  Here  t 


is  used  instead  of  t 
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alone  in  order  to  avoid  any  (distributional)  singularity  o£  tl>e  Green's 
function  at  the  integral  limit.  i'urtlior,  we  have  replaced  x”  and  x, 

■1  4 

by  the  equivalents  t and  t 

o 

The  second  integral  in  ( 78)  demands  special  care  in  i>erforming 
the  integration  over  the  time  coordinate  x°  since  v j is  an  explicit 
function  of  time,  so  that  the  integration  over  the  spatial  coord  inate.s 
is  to  be  performed  before  that  over  the  time  coordinate.  Further,  the 
resulting  integrand  for  the  time  integration  need  not  be  continuous, 
as  was  noted  earlier.  The  integral  in  question  is 


where  we  let  dt'^  ::  dx°  dx®  dx“  denote  the  purely  spatial  volume 

I element.  Now  we  observe  that,  for  a spatial  integral  of  a functioii 

! F of  the  medium  deformation  or  flow,  evaluated  over  a volume  v.'.| 

I 

!'  with  a function  of  time 

f 

_d_ 
dt 

, o 

1 
t 

! which  fol  lov/s  from  the  transpoT't  theorem  of  equat  Ion  (1). 

1 
I 

Here  the  surface  integral  is  due  to  mover'ent  of  the  boundary 
with  velocity  U . 'lliis  form  is  valid  no  m.itler 


/, 


F du°  = 


/ 


1?- 

Dt 


of  + 


f,(t  ) 
1 o 


f 


FU 


I id  a 
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what  the  value  of  U , that  is  whether  the  boundar>  of  move.; 

with  the  material  or  not.  Hence  it  is  immediately  applicable  to  the 
integral  being  considered,  where  we  may  regard  the  integrand  in  brtickets 
as  corresponding  to  K , so  tliat  : 


where  the  first  term  on  ttie  right  hand  side  must  bo  considered  as  a 
St ielj es  integral  (Archambeau,  1972;  Minster,  1973). 

In  the  theory  of  elastic  contlnua  and  in  linear  continuum  theory 
the  second  integral  in  (80)  is  neglected  since  the  boundary  movement 
is  with  the  particle  velocity  and  the  integral  term  is  small  - a second 
order  effect  due  to  transport  which  is  negligible  for  solids.  However, 
when  the  boundary  is  not  simply  a material  boundary  and  may  in  fact  move 
with  a normal  velocity  U-n  much  larger  than  the  particle  velocity,  then 
this  term  cannot  be  neglected  since  it  may,  in  fact,  be  as  large  as  the 


traction  term  at  the  boundary.  We  are  dealing  with  precisely  the  case  in 
whicii  ji  may  be  large  along  that  part  of  the  boundary  of  corresponding 

to  Z and  so  the  term  will  be  retained  for  Integration  over  the  "phase 
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boundaries  of  the  medium,  since  there  lias  the  value  of  the  particle 

velocity.  For  consistency  with  previous  neglect  of  second  order  terms 
we  are  in  fact  forced  to  neglect  material  boundary  position  variations 
for  which  in  (80)  is  equal  to  the  particle  velocity.  In  expressing 

(80),  however,  we  will  retain  for  the  surface  integral  limit,  but 

consider  to  be  zero  on  material  boundaries,  which  in  effect  neglects 

this  transport  Integral  term  for  boundaries  other  than  I . In  addition, 
U is  of  course  a spatial  function  so  that  its  value  may  be  negligible 
on  parts  of  the  surface  E as  well.  Indeed,  for  spontaneous  failure 
processes  we  would  expect  it  to  be  of  the  order  of  the  particle  velocity 
over  much  of  the  failure  surface  except  near  the  expanding  ’’edges*'  of  the 
failure  zone, which  usually  has  a elongated  shape.  The  magnitude  of  U 
on  I is  of  course  determined  by  the  jump  conditions  in  (32),  so  that  U 
is  not  arbitrary  but  determined  by  the  radiation  field  Itself.  The  shape 
of  tne  failure  boundary  is  thus  determined  by  the  spatial  variation  of  U , 
so  that  some  parts  of  E may  move  rapidly  while  other  sections  may  simply 
move  with  the  particles  as  a simple  material  boundary. 

Now  substituting  (79)  and  (80)  into  (78)  yields: 


(81) 
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From  previous  definitions  of  and  and  from  the  definition  of 

the  elastic-inertial  tensor  we  have 


"oB  "b 


r *1 

L o J 

[g™  - ^ iT*l 

[ kJi,  ^ at  '^s 


However  as  previously  noted,  we  must  neglect  terms  involving  products 
of  the  particle  velocities  in  comparison  with  traction  terms  or  tf>rms 
Involving  products  of  the  particle  velocity  and  • Therefore  we 

* 

can  replace  appearing  in  these  identities  by  alone.  Thus  in 

this  linearized  theory  we  use 


o 


ry-»ni  k 

^aB  "^B  " ^ 


and  observe  that  these  factors  appear  in  (81).  Further,  we  observe  that,  from  (74) 


r 

L 


9u 

.m  k 

"k  3t 


Using  these  identities  in  the  integrals  in  (81),  we  get 


t 

f — [g^  T - U 1 d^x  = f dt  f n.  da°  - 

J [^aoB  aaBj  o J o J B B 


a B 
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Thus  the  basic  representation  theorem  for  u^(X)  > expressed  by 


equation  (76),  is 


Ait  u^(X) 


-f 


pf  d^x 

an  o 


t'"  t-^  P 

^ f ^ L.  1 1 


dv 


9v 


o ^ u. 


(8A) 


Here 


tP  pP 

J„  = G T„-u  0 


3 a a3  ct  ap 


and 


= (n^.n^.n^,-  U"  n) 


Written  out  explicitly  in  three  dimensional  spatial  components,  this 
result  takes  the  form: 

.+ 


Attu^(X) 


c /^m  - 3 
p f,  G d X 
k k o 


o IT, 


/ L S 


3u, 


o / 


/ 3G“  3G™  \ 

“ “k  y^ijk«.  ^ 
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,m 


3u  n 
-,m  k 


k 3t 


dP 


(85  i 


o P, 


Here  p and  C..,  „ are  functions  of  the  spatial  coordinates  x, 
ijkx.  k 


alone. 
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6.  Approximate  Radiation  Field  Solutions  for  Growing  Failure  Zones 


,6 


The  fields  u and  G satisfy  the  boundary  conditions 
Y Y ^ 


B u = b ; xeSn 

ay  y a ' - 


B = 0 

ay  y 


XCd^ 


(86a) 


or,  equivalently 


.(1). 


"cxB  ^6 


= ^ 

aB 


xe8Q 


'’s  • ° 


(86b) 


; xe8(1 


where  the  superscripts  (1)  and  (2)  refer  to  the  quantities  evaluated 
by  limits  on  the  surface  approached  from  the  inside  or  from  the 
outside  of  n , respectively.  Thus  the  result  in  (84)  reduces  to: 

^ + 

/*  1 1 A i 

Attu 


(x)  = fp  f d^x  - f dt  f g’^[ 
— J a a o J o J a 


''e> 


Q 


/if 

n I-  If 


tM  . o 

J,  du 


(87) 
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(2) 

with  "prescribed"  on  the  boundary  of  . The  reduction 


of  (84)  to  (87)  is  a consequence  of  the  choice  of  the  Green's 
function,  which  eliminates  the  integral  involving  the  unknown 


displacement  field  on  the  boundary  of  , as  well  as 


from  the  fact  that  the  "generalized  tractions"  rig  appear 


explicitly  in  the  boundary  integral  over  so  that  the 


"jump  condition 


" fvs  "el  * " 


can  be  used  to  introduce  the 


action  of  the  material  external  to  at  the  boundary  of 


'1  • 


It  is  clear,  however,  that  determination  of  the  Green's  function 
satisfying  the  boundary  conditions  in  (86)  is  at  best  difficult 
and  quite  impractical  for  applications  in  view  of  the  necessity 
of  satisfying  boundary  conditions  on  the  growing  failure  boundary 


as  well  as  on  the  material  boundaries  of  the  medium.  We  can. 


,3 


instead,  introduce  a Green's  function  defined  on  a space 

a 


n'  that  includes  the  problc;7>  space  (2  , i.e.,  . In 

particular  we  can  take  Q'  to  be  the  space  bounded  by  the  same 


external  boundaries  as  Q but  without  the  internal  boundary 


,8 


created  by  the  failure  zone.  Then  is  given  by 


L r^(x  ; X ) 

ay  y — —o 


A®(x 

a — 


X ) 


xen' 


(88) 


X ) = 0 


xg9(2' 
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where  is  defined  throughout  as  if  the  whole  medium  were  in 

its  initial,  unfailed  state.  The  Green's  function  defined  by  (88)  has 

O 

the  same  properties  as  and  is  also  such  that  inner  products  over 

the  subregion  (2  of  Q'  can  be  defined.  Thus  the  Green's  representation 
of  (84)  is  reproduced,  with  G^  replaced  by  . However  the  price 

paid  for  the  presence  of  the  simpler  Green's  function  is  that  (84)  cannot 
be  reduced  to  the  simple  form  of  the  solution  given  in  (87).  In  particular, 
we  get,  from  (84)  and  (88) 


47r 


u (x)  = /*p  f d^x  - dt  f r^[ 

\i—Jaa  o J o J a 

n o 9u, 


/ •[/': 


du 


(89) 


where 

is 

the  difference  between 

the  spatial 

boundaries  of  (2 

and 

(2'  . 

With 

defined  as  above  this 

is  just  the 

failure  surface 

Y.  . 

Here 

^aB 

is 

the  "generalized  Green' 

s stress  tensor"  associated 

with 

r^‘ 

a 

, that 

is 

a o 

dx^ 


and  J,  is  obtained  from 
A 


J'*  = 


- u 

a aR  a aB 


1 


Thus  we  get  an  extra  integral  in  (89)  involving  u^^  on  the 
failure  boundary  E . Clearly  this  term  arises  from  the  fact 
that  we  have  in  effect  relaxed  the  boundary  conditions  required 
for  the  Green's  function. 

The  natural  boundary  conditions  on  Z do  not  involve  u^ 
directly,  but  instead  specify  conditions  relating  space  and  time 


derivatives  of  u across  I , i.e 
a 


’ _Ja6  %_,T 


0 . Hence 


u^  on  E is  unknown  and  (95)  is  an  integral  equation  for 
rather  than  a solution. 

Note  that  instead  of  the  Green’s  function  , we  can  use  the 

Y 

infinite  space  Green's  function  PjJ  . The  latter  affords  a relatively 
simple  closed  form  expression  (e.g.,  Maruyama,  1963).  In  that  case, 

I 

however,  the  boundary  integral  over  0 in  the  representation 

theorem  (89)  extends  over  all  boundaries  of  , including  exterior 

boundaries  as  well  as  E . In  addition,  we  must  consider  the  process 
to  occur  in  an  isotropic  and  homogeneous  space  from  the  onset.  This 
is  quite  an  acceptable  assumption,  to  first  order  at  least,  in  the 
vicinity  of  the  failure  zone. 

Successive  approximations  to  the  solution  of  (89)  can  be  obtained 
by  iteration.  In  this  particular  case  we  take  the  first  iterate,  > 

to  be  the  final  term  in  (89),  the  "relaxation"  , or  "initial  value"  term. 


Hence,  using  the  infinite  space  Green's  function  we  define: 


j o 


jy  dv° 

4 
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A correction  to  this  term,  which  yields  the  next  iterate  (x)  is 

defined  by 


4tt 


/•  r 

r “1 

uj(x)^  = 

■f  ‘"c,/”, 

r' 

^ ^ n 
a [_^a3  ^3 

da 


da 


(91) 


where  the  integration  is  only  over  the  failure  surface  I . (Generally 
the  volume  integral  over  the  body  force  density  f^  is  neglected  as 
being  small  and  is  not  included.)  In  the  past  we  have  referred  to  ii^(2t) 
as  the  "transparent  source  approximation  to  the  direct  source  field", 
since  it  neglects  the  scattering  at  the  failure  boundary.  The  field 
jj^^(3t)  will  simply  be  referred  to  as  the  "approximate  direct  source 
field". 

By  adding  a general  eigenfunction  expansion  with  adjustable 
coefficients,  which  is  regular  every\>;here  in  the  domain  Q (or  Q'  ); 
one  can  then  satisfy  boundary  conditions  on  the  other  boundaries.  For 
example,  for  a problem  involving  failure  in  the  Earth,  we  could  use  a 
layered  spherical  model,  and  hence  a vector  spherical  harmonic  expansion 
for  each  layer  zone.  One  could  also  express  the  Green's  function  in 

(89)  as  an  eigenfunction  expansion  and  use  it  in  place  of  F^^  in  (90) 
and  (91)  to  obtain, directly, the  approximate  source  fields  in  the  domain 


Q'  . 


r 


It  is  clear  that  the  relaxation  or  initiai  value  term  plays  a 
dominant  role  in  this  procedure.  In  fact,  we  have  pointed  out  that  it 
represents  the  essentials  of  the  elastic  wave  radiation  attributable  to 
the  source.  Therefore  we  will  consider  its  expression  in  more  explicit 
and  simplified  form  for  spontaneous  processes  such  as  failure  under  ck, 
static  or  quasi-static  initial  stress  condition.  Related  investigations  have 
also  been  considered  by  Archambeau  (1968,  1972)  and  Minster  (1973). 

If  we  consider  the  field  u^(2{)  at  any  particular  time  tj^  we 
can  represent  its  evolution  at  any  later  time  t in  terms  of  the  initial 
value  integral  (e.g..  Love,  1944) 


4tt  u (x,  t;  t ) 
m — 1 


Ui(ti) 


3u 

..m  k 


o 1 


This  expression  in  the  radiation  field  due  to  all  effects  occurring  up 
to  time  tj^  , including  prior  rupture  growth.  It  is  therefore  the 
rgdiation  that  would  be  observed  for  t > tj^  if  the  failure  xone 
boundary  were  fixed  for  times  ^ ^ consider  the  same 


initial  value 


representation  of  the  field  at  t^  + 6tj^  , where  some 


additional  change  in  the  failure  boundary  has  occurred  in  the  time 
increment  6tj^  , then  the  difference  in  the  radiation  fields  is  given  by 


u^(x,  t;  t^  + 6t  ) - u^(x,  t;  t,)  = 
m — 1 1 m — 1 


u^(tj^  + 6t^r 


3t  / 3t 
. o / o 


dv°  + 0(6tJ) 


t = t, 
o 1 
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where  -5--  6t.  is  the  incremental  initial  value,  equal  to  the 

i 

o t =t, 
o I 

change  in  the  equilibrium  field  caused  by  incremental  growth  of  the  failure 
zone,  and  where  we  have  used  the  fact  that  accelerations  remain  finite 
so  that  the  velocity  is  continuous. 

Dividing  through  by  6tj^  , and  taking  the  limit  as  6tj^  approaches 
zero  gives 


^1) 

3t, 


vpti)' 


3t  / 9t 


t = t, 
o 1 


For  failure  occurring  at  a finite  rate  and  of  total  duration  T , then 
the  derivative  with  respect  to  t^  vanishes  for  t^  > T . Thus 
integration  of  (92)  yields 


'>  ■ c(i.  t)  + / dt  J 

J ^ \ ^ 0/  o 


t = t, 
o 1 


where  the  constant  of  integration  is  fixed  by  the  choice  of  reference 
equilibrium  state.  If  we  chose  to  measure  u^  relative  to  the  final 
equilibrium  state  after  the  complete  boundary  has  been  formed,  then 


it 

c(x,  t)  = t)  - '1^) 


so  that 


it  it 

u (x,  0)  = u (x,0)  - u (x,  t) 
m — m — m ~ 
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If  wc  refer  tlie  d isplacements  to  Lliu  iiiiLial  equilibrium  state  tiien 

c(>c,  t)  = 0 . Iti  either  case  c(x,  t)  vanishes  for  t > l . i li  i s 

result  is  formally  exact;  if  or  f'^  is  used,  however,  source  boundary 

a a 

scattering  is  uegiecteJ.  On  the  other  hand,  the  equilibrium  field 
changes  expressed  by  t^)  account  for  both  the  spatial  and 

temporal  behavior  of  the  failure  boundary,  so  that  the  only  term 
neglected,  as  expressed  approximately  in  (91)  represents  the  dynamic  inter- 
action of  the  radiation  field  and  the  boundary,  that  is  scattering. 

We  expect  this  scattering  to  be  significant  only  for  wave  lengths  smaller 
than  the  characteristic  source  dimension.  This  is  confirmed  by  the 
results  of  Burridge  (1975),  Koyama  ot  al . (1973)  and  Minster  and  Suteau 
(1976).  These  authors  considered  expanding  spherical  failure  zones 
in  an  uniformly  prestress  infinite  space.  Scattering  by  E is  shown 
to  affect  the  spectrum  only  for  frequencies  greater  that  the  source 
characteristic  frequency  f = b'  /L  , where  L is  the  maximum  source 

dimension  and  U the  average  failure  ratt;  ind  the  effects  are  not 
K 

large.  In  the  time  domain,  the  "transparency"  approximation  leads  to  a 
minor  shortening  of  the  far-field  pulse  and  elimination  of  its  small 
amplitude  oscillatory  coda. 

7 • Conclusions  and  Consequences 

Some  of  the  conclusions  that  we  have  dra\m  from  the  results 
of  the  work  described  are: 

(i)  Conservation  relations  expressed  on  an  expanding  material 
transition  surface  (or  failure  zone  boundary)  show  that  traction  and 
particle  velocity  Jumps  across  such  a boundary  cannot  be  specified 
independently,  as  has  been  (implicitly)  assumed,  for  example,  in 
dislocation  and  stress  pulse  (crack  theory)  models  for  failure.  In 
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particular,  it  has  been  shown  that  the  proportionality  factor  between 

particle  velocity  clianges  and  traction  cltanp.es  across  tiie  failure 

boundary  is  pU  , wltti  U tlie  rupture  rate  and  p the  initial 
R R 

material  density. 

(2)  Introduction  of  a transition  energy  function  provides  a 
metins  of  characterizing  failure  processes  in  terms  of  energy  requirements. 
Tills  function  can  be  determined  from  seismically  estimated  rupture  var- 
iables— in  particular  from  the  stress  drop  (on  the  failure  boundary) 

and  the  rupture  velocity.  Hence  the  transition  energy  (density) 
function  can  be  obtained  for  natural  failure  processes  and  compared 
to  laboratory  determined  values  which  sliould  lead  to  a more  complete 
quantitative  description  and  subsequent  understanding  of  failure  processes 
in  the  eartli.  (See  also  Hussein!  et  al.,  1975.) 

(3)  Since,  rupture  rates  have  been  shown  to  be  directly  proportional 
to  stress  drops  on  the  failure  boundary  and  also  proportional  to  the 
change  in  the  material  moduli,  it  follows  that  observations  of  high 
rupture  rates  for  earthquakes,  near  the  shear  velocity  of  the  medium, 
imply  near  total  loss  of  shear  strength  and  probable  melting  during 
failure.  (See  also  Archambeau,  1976.) 

(4)  The  Green's  function  solution  for  the  radiation  field 
exterior  to  the  failure  zone  shows  that  the  field  depends  explicitly 
on  the  failure  boundary  growth  rate,  whereas  previous  radiation  field 
representations  did  not  contain  any  such  explicit  dependence  since 
dynamical  boundary  conditions  were  ignored  in  favor  cf  heuristic 
kinematical  representations  of  boundary  growth  effects.  This 

suggests  that  source  representations  such  as  those  generated  by  displace- 
ment dislocations  and  stress  pulse  equivalents,  that  are  formulated 
2ntirely  in  terms  of  equivalent  applied  boundary  displacements  or 
equivalent  applied  boundary  tractions,  must  be  carefully  reconsidered 
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since,  at  present,  they  do  not  account  for  the  boundary  integral  terns 
involving  the  surface  growth  rate  (i.e.,  the  surface  integral  terns 
with  integrands  involving  explicitly,  as  given  for  example  in 

equation  (85)  ). 

(5)  The  Green's  function  solution  shows  that  an  appropriate  first 
order  representation  of  the  radiation  field  is  obtained  using  the 
"initial  value”  or  "relaxation”  terra  in  the  Green's  integral  representation. 
The  radiation  field  predictions  based  on  this  integral  representation 
can  be  directly  related  to  the  physical  parameters  for  the  failure  process 
(i.e,,  the  transition  energy  and  the  material  rheology  moduli  before 
and  after  failure  as  well  as  the  rupture  rate  and  stress  level  changes) 
and  can  be  used  for  the  representation  of  the  radiation  from  complex 
failure  phenomena  in  inhomogeneous  media  that  is  also  inhomogeneously 
prestressed.  Since  strong  prestress  inhomogeneities  are  considered 
to  be  present  for  most  if  not  all  failure  processes  (especially  for 
earthquakes)  and  to  have  a very  important  effect  on  the 
character  of  the  radiated  elastic  wave  field,  then  it  follows  that  the 
ability  to  theoretically  predict  the  field  from  such  complex  sources 
is  critical  to  the  advancement  of  understanding  of  the  mechanics 


of  failure  processes,  particularly  in  the  earth. 
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